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ABSTRACT 



Most neighboring stars are still detected as point sources and are beyond the angular res- 
olution reach of current observatories. Methods to improve our understanding of stars at 
high angular resolution are investigated. Air Cherenkov telescopes (ACTs), primarily used for 
Gamma-ray astronomy, enable us to increase our understanding of the circumstellar environment 
of a particular system. When used as optical intensity interferometers, future ACT arrays will 
allow us to detect stars as extended objects and image their surfaces at high angular resolution. 

ACTs are used in gamma-ray astronomy to investigate violent phenomena in the universe. 
However, this technique can also be used for stellar astrophysics on some isolated sources. Such is 
the case with the X-ray binary LS I -|-61°303 which was detected in the TeV range. A gamma-ray 
attenuation model is developed and applied to this system. This models allows us to place 
constraints on fundamental properties of the system. However, a much better understanding of 
this system, and more so of nearby bright stellar systems, could be obtained with high angular 
resolution techniques. 

Optical stellar intensity interferometry (SH) with ACT arrays, composed of nearly 100 
telescopes, will provide means to measure fundamental stellar parameters and also open the 
possibility of model-independent imaging. A data analysis algorithm is developed and permits 
the reconstruction of high angular resolution images from simulated SH data. The capabilities 

and limitations of future ACT arrays used for high angular resolution imaging are investigated 
via Monte-Carlo simulations. Simple stellar objects as well as stellar surfaces with localized hot 
or cool regions can be accurately imaged. 

Finally, experimental efforts to measure intensity correlations arc expounded. The function- 
ality of analog and digital correlators is demonstrated. Intensity correlations have been measured 
for a simulated star emitting pseudo-thermal light, resulting in angular diameter measurements. 
The StarBase observatory, consisting of a pair of 3 m telescopes separated by 23 m, is described. 
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CHAPTER 1 

HIGH ANGULAR RESOLUTION 
ASTRONOMY AND COHERENT 
LIGHT 

1.1 Introduction 

Most stars are still detected as unresolved point sources at optical wavelengths, and our 
understanding at high angular resolution results from physical modeling of their light spectrum 
and variability. Stellar surface brightness distributions can be predicted to some degree, but 
stellar atmospheres (may) have convection zones, mechanically driven matter flows, radiation 
driven winds, accretion, and other complex phenomena that are difficult to investigate without 
the recourse to high angular resolution imaging. Imaging surface structures at near-optical 
wavelengths may provide direct evidence for many of these effects and is another mean to test 
our current understanding of stellar atmospheres and stellar evolution. However, it is only for a 
few stars that actTial images have been obtained using different techniques. 

Recent results of optical high angular resolution astronomy, particularly Michelson interfer- 
ometry, have started to reveal stars as extended objects and have increased our understanding of 
effects such as those mentioned above. Reconstructed images of stars with non-uniform radiance 
distributions have, in some cases, confirmed our understanding of stars, but have also surprised us 
in others, reminding us that we have still much to learn. However, most stars are still beyond the 
angular resolution of current methods, which is why we propose to revive an extremely successful, 
yet abandoned technique, namely Stellar Intensity Interferometry (SH). Being insensitive to 
atmospheric turbulence and imperfections in optics, a previous SH experiment pioneered by 
Robert Hanbury Brown (Hanbury Brown, 1974) produced more scientific results in less than a 
decade than several amplitude (Michelson) experiments combined several decades later. 

Over the course of this work, we have realized that the techniques used to study stars are 
just as interesting as stars themselves, which is why some historical and technical background 
is given in Chapters 1-3. SH will be compared to other existing techniques such as Michelson 
interferometry and speckle interferometry, which allows one to obtain data that is similar to 
intensity interferometry (Fourier magnitude data). We will discuss how interferometric data 
allows one to obtain high angular resolution images, and how these techniques typically allow 
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to measure the Fourier transform of the radiance distribution of the star. The main challenge 
for recovering images of stars with SIX is the lack of phase information, which is discussed in 
Chapter 4. Some new results in phase retrieval techniques were obtained with a method relying 
on the theory of analytic functions. Provided sufficient Fourier coverage and signal-to-noise, this 
method will allow to obtain high angular resolution images of stars. 

On the other side of the photon energy scale, ground-based gamma-ray astronomy, discussed 
in Chapter 5, is a flourishing field which has allowed to detect very high energy (TeV) radiation 
emitted from galactic and extragalactic objects. In Chapter 5, a short detour from high angular 
resolution astronomy is taken to analyze a particular high energy emitting binary system {LSI+ 
61° 303) consisting of a fast rotating main sequence star experiencing mass loss and a compact 
object whose nature is still subject of debate. This object displays a TeV light-curve that shows 
a modulation with the same periodicity as the binary, and almost begs for the development of 
a toy model in terms that describes the 7-ray attenuation. When fitted to the data, this model 
allows us to constrain some fundamental parameters of the system (Nunez et al., 2011). However, 
many of the constraints that can be placed with TeV observations are still debatable, and optical 
high angular resolution data may answer many questions on the nature of this object. In fact, 
the optical requirements for optical SII are adequately met with the large light collectors used in 
7-ray astronomy, and this has prompted a recent revival of optical SII using lACTs (Le Bohec 
&; Holder, 2006). 

The recent success of 7-ray astronomy has motivated the construction of large lACTs. In 
the case of the Cherenkov Telescope Array project (CTA), it is anticipated that there will be 
nearly 100 telescopes that will be separated by up to 1km (Consortium, 2010). When used 
as an intensity interferometers, CTA will complement the science done with existing amplitude 
interferometers by observing at shorter wavelengths 400 nm) and increasing the angular 
resolution by nearly an order of magnitude. In Chapters 6-7, the capabilities of lACTs used for 
SII will be quantified via simulations, and the most important result of this work was to prove 
that high angular resolution images can be obtained from optical SII data obtained from these 
arrays by applying the phase retrieval techniques discussed in Chapter 4 (Nunez et al., 2012a). 

There are also ongoing efforts for performing an intensity interferometry measurement. The 
final chapter (8) starts with a discussion of our current experimental efforts to measure intensity 
correlations from thermal sources in the laboratory and from stars with a pair of small Cherenkov 
telescopes. The detection of intensity correlations from a thermal source has proven to be elusive 
due to the extremely short coherence times and shot-noise dominated data. In order to gain some 
experimental understanding, intensity interferometry data were obtained by using a light source 
with a much longer coherence time by breaking the coherence of laser light (pseudo-thermal 
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light). Here, angular diameters of small pinholes emitting pseudo-thermal light were measured. 

1.2 A brief history of high angular resolution 
astronomy and intensity interferometry 

The problem of resolving stars dates back at least 400 years. Even before considering the 
problem of imaging, a very old problem is that of measuring the angular diameter of stars. 
Galileo placed an upper limit on the angular size of Vega by measuring the distance he needed 
to be from a string in order for the star to be obscured (Galilei, 1953). Galileo's reported angular 
size of 5 arcseconds is most likely a measure of atmospheric scintillation. Newton also estimated 
the angular size of stars by assuming they were like distant suns and found the angular diameter 
of a first magnitude star like Vega to be 2 milliarcseconds (mas), which is actually very close to 
the accepted value of 3 mas. 

By the late 1800's it was already known that the maximum angular resolution was not 
ultimately determined by the size of the aperture, but by the atmospheric turbulence or "seeing" 
of ~ 1 arcsecond. This limits the practical size of a telescope to about 10 cm for resolving small 
objects. Interferometry was what really revolutionized high angular resolution astronomy when 

H. Fizeau (1868) proposed to mount a mask on a telescope to observe interference fringes. This 
was implemented by M. Stephan in the Marseilles observatory with an 80 cm telescope, the 
largest instrument at the time (1874). Since they clearly saw interference fringes on every star 
they pointed at, they provided an upper limit of 0.16 arcseconds for the angular diameters of 
stars^ and improved the angular resolution by an order of magnitude. The angular resolution 
would still have to be improved by several orders of magnitude in order to detect stars as 
extended objects. The first successful measurement of angular diameters was performed by 
Michelson and Pease (1920) at Mt. Wilson, and this is considered the birth of high angular 
resolution astronomy. In his preliminary experiments with a 36" telescope, Michelson measured 
the angular diameter of Jupiter's Galilean moons (Michelson, 1891). With an interferometer 
with a maximum baseline^ of 6 m, they first measured the angular diameter of the red giant 
Betelgeuse to be 47 mas (Michelson & Pease, 1921), and in total measured the angular diameter 
of 6 giant stars, all within tens of milliarcseconds. 

Michelson was aware of one of the main difficulties in optical stellar interferometry, namely 
that of atmospheric turbulence effects. In his preliminary experiments, he noticed that fringes 
would drift (jitter) in time as a result of turbulence induced path differences between the two 

^The contrast of fringes is smaller for sources that subtend a larger angle in the sky as will be shown in section 

I. 3.1. 

^We will use the term "baseline" to refer to the separation between points where the light signal is received. 
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interfering beams. Fringes can drift on timescales of the order of milliseconds, and Michelson's 

work was truly remarkable considering that he measured the fringe contrast (visibility) with 
the unaided eye. Pease (1930) then built a larger 15 m version of their previous interferometer, 
but this project failed to give reliable results and was later abandoned (1937). This "failed" 
experiment clearly shows another main difficulty in Michelson optical interferometry, namely 
the need for very high (subwavelength) precision optics. It took over three decades for the next 
breakthrough to occur in optical interferometry. 

The understanding of interference phenomena in terms of coherence theory was a next crucial 
step for the advancement of high angular resolution astronomy. Michelson never mentioned the 
word "coherence", but he knew how the fringe visibilities should behave as a function of the 
baseline for uniform disks, limb-darkened stars and binaries. The first investigations of coherence 
phenomena are due to Verdet (1865) (Mandel & Wolf, 1995), when he used the angular diameter 
of the sun to estimate the (transverse) coherence length of the sun to be 0.1mm (see Section 
1.3.1). Between the 1920s and the 1950s several notorious scientists, including Weiner, van 
Cittert, Zernike, Hopkins and Wolf, participated in further investigations of coherence theory. 
The outcome of these studies was a method to quantify the correlations between fluctuating fields 
at two space-time points and the dynamical laws which correlations obey in free space (Mandel 
& Wolf, 1995). A particularly important result was to relate the degree of coherence, partially 
obtained from the visibility of interference fringes, to the Fourier transform of the radiance 
distribution of the star. This last statement is known as the van Cittert-Zernickc theorem 
(section 1.3.4), and led to the development of aperture synthesis, which was first applied in 
radio astronomy by M. Ryle (1952). Radio interferometry became a flourishing field and is still 
responsible for most of the highest angular resolution images available today. The techniques 
learned during this period along with improved understanding of coherence theory quickly led 
to the development of intensity interferometry. 

Intensity interferometry was born in 1949 as a result of attempts to design a radio inter- 
ferometer that could measure the angular scale of two of the most prominent radio sources, 
Cygnus A and Cassiopeia A. R. Hanbury Brown conjectured that the electronic noise recorded 
at two different stations was correlated, that is, that low frequency intensity fluctuations were 
correlated. At the time, Hanbury Brown thought that the main advantage of this technique 
over conventional radio interferometry was that it did not require synchronized oscillators at 
two distant receivers so that observations could be done with at longer baselines. With the help 
of R. Twiss, a formal theory of intensity interferometry was developed (Brown &; Twiss, 1957) 
(Chapter 3), and before long, they built the flrst radio intensity interferometer and eventually 
succeeded to measure the expected correlations and the angular sizes of Cygnus A and Cassiopeia 
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A (1952). They realized however that they had "used a sledge-hammer to crack a nut" since the 
baseline was only of a few kilometers, and this could have been accomplished with conventional 
radio intcrfcronictry. Nevertheless, this experiment made them realize one of the most important 
advantages of intensity interferomctry: they noticed that even when the sources were scintillating 
violently due to poor atmospheric conditions, the correlation was not affected. Correlations of 
intensity are essentially unaffected by atmospheric turbulence (see Section 6.1 for more details). 

In order to test the theory at optical wavelengths, R. Hanbury Brown and R. Twiss performed 
a series of laboratory experiments to measure the expected correlations (Brown & Twiss, 1958). 
In these experiments, two photomultiplier tubes were placed a couple meters from a pinhole 
that was illuminated by a narrow-band thermal source, and the correlation was then measured 
as a function of detector separation. These experiments were surrounded by controversy, and 
were followed by several experiments which seemed to contradict Hanbury Brown's and Twiss's 
results (e.g., Adam et al. (1955), Brannen & Ferguson (1956)). Most of these experiments used 
broadband light sources, which correspond to low spectral densities, and would have needed 
extremely long integration times in order to detect significant correlations (Hanbury Brown, 
1974). These experiments were confirmed several times (e.g., Morgan & Mandel (1966)), and the 
subsequent quantum understanding of Hanbury Brown's and Twiss's work led to the development 
of quantum optics. 

The following step was to construct an optical stellar intensity interferometer. Since large 
light collectors were needed, Hanbury Brown was concerned about the cost of two large optical 
telescopes. He soon realized that, for the same reason that measurements were insensitive to 
atmosphere induced path delays, there was no need for high precision optics, and "light bucket" 
type detectors could be used. In the small town of Narrabri (Australia), two movable reflectors 
were placed on circular tracks (see Figure 1.1) with a control station at the center were signals 
were correlated. The circular track allowed the baseline to be maintained perpendicular to the 
distant star. Between 1965 and 1972, the Narrabri stellar intensity interferometer measured the 
angular diameter of 32 bright stars, a considerable extension to the long standing catalog of 6 
stars obtained previously by Albert Michelson. 

The main work at Narrabri was completed when the diameters of the brightest stars in the 
southern hemisphere were measured, and the sensitivity limit of the instrument had been reached. 
Plans for a larger instrument were abandoned when it was shown that the same sensitivity could 
be reached with amplitude interferometry in 1/40 of the time (Labeyrie et al., 2006). Even so, it 
took over 30 years for amplitude interferometrists to overcome all the difficulties and to match 
the sensitivity of intensity interferometry. Even though stellar intensity interferometry was not 
further pursued, the same technique has been applied in other fields such as high energy particle 
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Figure 1.1. Light collectors used in the Narrabri stellar intensity interferometer. Image from 
Sky & Telescope, vol. 28, pp. 2-7, 1964. 

physics to probe nuclei at very small scales^ (Baym, 1998). 

An interesting technique developed during this time by A. Labeyrie (1970) was that of 
speckle interferometry (Labeyrie, 1970), where short exposures of a single telescope reveal 
atmospherically induced speckle patterns. Speckle patterns can be thought of as being generated 
by a random mask that changes at short (ms) time-scales, where each subaperture has a typical 
size (~ 10 cm) given by the atmospheric seeing. The interpretation of these speckle patterns 
allowed to obtain diffraction limited information as opposed to being limited by atmospheric 
seeing. This technique was used to obtain results on hundreds of fainter sources than those 
previously observed by Michclson, many of which were found to be binary stars. ^ 

Labeyrie's work started a revival in optical stellar (amplitude) interferometry which resulted 
in the first successful beam combination of two telescopes separated by 12 m (Labeyrie, 1975). In 
this case, Labeyrie saw fringes on the bright star Vega, and soon after, several projects developed 
the field that is now known as optical long baseline interferometry. AmplitTidc interferometry 
is currently producing very interesting results since technological advancements have made it 
possible to have baselines of the order of ^ 100 m and to combine beams from more than two 



^In this case by measuring correlations of particles with integer spin such as pions. 

*Many of these observations were carried out by CHARA (Center for High Angular Resolution Astronomy) 
between 1977 and 1998, and the success of this program laid the groundwork for CHARA's entry into long baseline 
optical interferometry (Hartkopf & Mason, 2009) 
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telescopes. This has resulted in several very high angular resolution images of rapidly rotating 
stars such as Vega (Aufdenberg et al., 2006) and Altair (Monnier et al., 2007), and also the 
interacting binary /3-lyrae (Zhao et al., 2008), and the eclipsing binary e-aurigae (Kloppenborg 
et al., 2010). 



Interference phenomena are at the heart of high angular resolution astronomy observation, 
and therefore we start with a discussion on interference and its interpretation in terms of 
correlation theory. A plane wave can be formed from a point source located at infinity, and the 
study of interference phenomena allows us to measure deviations from a source being point-like. 
Light which has the same frequency and phase in all space is said to exhibit spatial coherence, 
and light whose frequency and phase are well defined^ over time is said to exhibit temporal 
coherence (Mandel Sz Wolf, 1995; Labeyrie et al., 2006). Interference can be clearly seen with 
plane monochromatic waves, which maximally exhibit spatial and temporal coherence, and the 
extent to which these phenomena can be observed is determined by the departure from these 
ideal conditions. 

1.3.1 A classical discussion on coherence 

The amplitude of the electric field E{t) of light at a fixed position can be expressed as the 
superposition of many monochromatic plane waves as 



If the electric field is non-zero during a time interval At, then E{u;) (which may be complex) 
has to have variations over a bandwidth Aa; comparable to variations of e*'^^* (Landau &; Lifshitz, 
1975), which occur when Aa; ~ 1/At. That is 



The time scale At is known as the coherence time, and is the time during which light can be 
considered approximately monochromatic with a defined phase. The length associated to this 
time is known as the coherence length given by 



1.3 Theory of partially coherent light 




(1.1) 



AuAt ~ 1. 



(1.2) 



well defined frequency corresponds to a field expressed as a single plane wave as opposed to a superposition 
of plane waves of different frequencies. 
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A similar argument can be made by expressing the electric field E{x) at a particular time as 
a superposition of plane waves emitted from different angular positions ^ of a far away source, 
i.e., 

E{x) = J E{e)e'^''^de. (1.4) 

If the electric field is different from zero over a distance Ax, then E{6) has to have variations 
over angles A^ comparable to variations of e^^^^^ , which occur when A0 ~ l/fcAx. That is 

A;AxA^~l. (1.5) 

The length Ax is known as the transverse coherence length and is the length in which the 
wave is approximately planar. 

Another quantity of interest is the coherence volume, given by 

Ay-Ax^A/-^. (1.6) 

The physical significance of this volume will now be explained. 

1.3.2 Quantum interpretation of coherence 

It is interesting to note the connection with quantum mechanics by finding the volume of 
phase space AV given by the uncertainty principle, that is 

~ — (1-7) 

Ap^.Apj,Ap^ 

Within this volume, photons arc indistinguishable, i.e. they belong to the same mode. If a 
source has a (small) angular size A0, then Ap^, = pO j6, and Apz ^ Therefore, 

This is equal to the coherence volume (eq. 1.6) found with classical arguments. In the 
quantum interpretation, the position of a photon is not defined within the coherence volume 
before a measurement takes place, and therefore interference phenomena can be seen. The 
classical example is given by the double slit experiment, where the wave-function can be described 
by the superposition of the photon having one position, or the other position, and the squared 
modulus yields a probability distribution with interference terms. In the quantum interpretation, 
the coherence volume is the region in which photons are indistinguishable (Mandel &: Wolf, 1995). 
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Even though they are indistinguishable, photons within the coherence volume are correlated (and 
not independent), since only symmetric states can occur between them. Brown and Twiss (1957) 
summarized the above by realizing that the physical significance of the coherence volume is that 
since photons within this volume can interfere®, they are indistinguishable. 

1.3.3 Orders of magnitude for transverse coherence 
lengths and coherence times 

As a first example, we can consider the sun, which is obviously not a point source, so we 
should not expect to detect spatially coherent plane waves. The angular size of the sun is ~ 0.5°, 
so the transverse coherence length of the sun at A = 500 nm is 0.05 mm. This means we could 
observe interference fringes from the sun if slits were separated by less than 0.05 mm. The most 
giant stars have angular diameters of tens of milliarcseconds, so that their transverse coherence 
length is of the order of tens of meters at A = 500 nm. Most stars have angular diameters of the 
order of ~ 0.1 mas, so that their transverse coherence length is of the order of a few hundred 
meters, hence the need for long baseline optical interferometry. 

The time resolution of fast photo-multiplier tubes is of the order of ~ 10~^s, which is several 
(3-4) orders of magnitude longer than the typical coherence time of a thermal light source, and 
poses a challenge for measuring temporal coherence effects (see Section 3.1). Very coherent light 
sources such as lasers on the other hand, have much longer coherence times of the order of 10~^ s. 

1.3.4 Interpretation in terms of correlations 

More understanding of interference phenomena can be gained by studying the correlation 
between partial beams of light that are separated in space and time. If two partial beams are 
superposed with a time delay between them that is much less than the coherence time, then the 
two are highly correlated. When a time delay of the order of the coherence time is introduced, 
the two beams no longer have the same frequency and phase, and therefore there is no correlation 
(Mandel & Wolf, 1995; Labeyrie et al., 2006). 

A similar reasoning can be made with spatial coherence. A distant extended object, such as a 
star, can be thought of as a collection of points emitting spherical plane waves. The light emitted 
from different points in the object is statistically independent and not correlated. However, if 
light from the object is being detected at two points in space, the superposition of the light 
emitted from all points of the source is more correlated when the two detection points are close 
together. That is, two detectors that are close together receive essentially the same light signal. 



It is important to emphasize that photons do not interfere with each other, but with themselves, due to the 
probabiUstic nature of their quantum state. 
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These statements can be made more quantitative by introducing the complex degree of coherence 
between the electric field of light E{f,t) at space-time points (ri,ti) and (r2,t2) as 

7(ri,ti;r2,i2) = — , (1.9) 
^{\E{r,M)?)mr2M)?) 

The brackets (...) indicate time averaging. The numerator in this expression arises naturally 
as an interference term when calculating the intensity of the superposition of two beams of light. 
The denominator is for normalization purposes. The complex degree of coherence is of central 
importance in interferometry, since it can be related to measurable quantities. In the case of 
amplitude interferometry, the complex degree of coherence can be related to the fringe visibility 
which is defined as 

Imax ~ Imin q 

T A- T ■ \ ■ ) 

■'■max ~ ^mm 

In fact, if two monochromatic waves of the same frequency and with amplitudes Ai and A2 
are combined, it is straightforward to show that the fringe visibility is^ (Labeyrie et al., 2006) 

^ = ^42^l7(n,i;r2,t)|. (l.ll) 

The visibility is the main observable in amplitude interferometry and allows us to measure 
the light coherence. Moreover, a measurement of the complex degree of coherence allows one to 
obtain information about the angular radiance distribution of the source as we will now show. 

1.3.5 Connection between source structure 
and light coherence 

We will consider a monochromatic light source and choose to set the origin of coordinates at 

the light source. Points on the source are labeled by positions x', and the position of a point on 

the source with respect to a far away observer is 

f=x-x'. (1.12) 

The observed electric field at x is 

E{x, t) = J ^^e'^'-cfx', (1.13) 

where k is the wave number 27r/A, and we have omitted the time variation e"^* and the random 
phases induced by turbulence among other factors. Now we make the approximation 



Eq. can be derived by first finding tlic intensity of tlie superposition of electric fields, finding the maximum 
and minimum intensities, and noting that cross terms in the electric field correspond to the degree of coherence. 
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\r\ f=i \x\ - x' ~ (1.14) 



so that 



E{x,t) = J A{x') exp likx' ■ — jdV. (1.15) 
Now we calculate the time averaged correlation {E{xi,ti)E*{xj,tj)) for the particular case 



when ti = tj . We should first note that 



{A{x')A{x")) = I{x')6{x' - x"), (1.16) 

where I{x') is the light intensity at point x'. This is because separate points on the source are 
not correlated over large distances. Now the time averaged correlation is 

{E{xi)E*{xj)) = C J /(f) exp (^f • ^ - f • I (fx', (1.17) 

where C is a constant. When ^ and <C \xj\, we can write the angle 9 as 

(1.18) 



x' x' 



We can now express the correlation as 



{E{xi)E*{xj)) = J I{e)e-'''^<^'-^^W^9 (1.19) 
= i{xi-Xj), (1.20) 

where I(zi — Zj) is the Fourier transform of the radiance distribution of the star. The Fourier 
transform goes from angular space to detector separation space. 

We have just shown that the complex degree of coherence is the normalized Fourier transform 
of the angular intensity distribution of the source (Labeyrie et al., 2006). This last statement 
is known as the Van Cittert-Zernike theorem. Therefore, by measuring the complex degree of 
coherence, or a related quantity, we can learn about the source structure. In Section 1.3.4 , 
we saw that the magnitude of the Fourier transform can be directly measured in amplitude 
interferometry with visibility measurements, and we will see that it can also be measured with 
intensity correlation measurements (Section 3.2). The phase of the Fourier transform is more 
elusive, and techniques have been developed for its measurement with amplitude interferometry 
(Section 4.3) and its recovery with intensity interferometry (Section 4.5). 



CHAPTER 2 



STELLAR ASTROPHYSICS AT HIGH 
ANGULAR RESOLUTION 

Most stars are still merely detected as point sources and not as the extended and complex 
objects they truly are. One can only speculate as to what will be learned once more is known 
of this high angular resolution world. An analogy can be made with extragalactic astronomy, 
where one could ask what would be the status of the field if galaxies were regarded as unresolved 
point sources. Some fundamental stellar parameters as well as several interesting effects can be 
studied with high angular resolution astronomy. In this chapter, a few representative topics of 
high angular resolution astrophysics are discussed. 

2.1 Angular diameters 

Stars can be characterized by measuring basic parameters such as the effective temperature, 
luminosity, chemical composition, mass and radius. A broad sample of these parameters provides 
effective constraints on stellar evolutionary models. Some of these quantities can be measured 
using traditional astronomical techniques. For example, the chemical composition can be mea- 
sured by studying spectral lines. The effective temperature can be obtained with knowledge 
of the integrated light flux and the physical radius, which can be obtained by measuring the 
angular radius and the parallax (distance). In order to constrain the position of stars in 
the Hertzprung-Russell diagram, radii measurements with uncertainties of a few percent are 
necessary (Aufdenberg et al., 2005). 

The angular extent of stars is typically less than 1 milliarcsecond, and is only tens of 
milliarcseconds for even the most nearby giant stars. In Figure 2.1, the estimated angular 
diameter for 35000 stars in the JMMC (Lafrasse et al., 2010) catalog are shown, and we can see 
that the high angular resolution world really starts below ~ 10 mas. Light received from most 
stars has transverse coherence lengths of several hundred meters (Section 1.3.3), so that model 
independent measurements of the physical size (at optical or near-optical wavelengths) can only 
be obtained through interferometric techniques (Labeyrie et al., 2006; Ten Brummelaar et al., 
2009), and there are still very few measurements of this basic fundamental parameter. At the 
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Figure 2.1. Number of stars for each estimated angular size from the JMMC stellar diameter 
catalog. The approximate angular resolution at (near) optical wavelengths for NSII, CHARA, 
and the Hubble instruments is shown for reference. 

time of this writing, 242 stellar diameters have been directly measured, out of which only 24 
correspond to main sequence stars (Boyajian et al., 2011). 

One of the main applications of angular diameter measurements is the understanding of stellar 
atmospheres. As was mentioned before, interferometric measurements can provide information 
of the effective temperature, and when combined with spectroscopic and spectrophotometric 
measurements, it can be used to create consistent models for stellar atmospheres (Labeyrie 
et al., 2006). The challenge for stellar atmosphere modelers is to reproduce stellar line spectra 
from knowledge of these fundamental stellar parameters. This approach uses an angular diameter 
measurement, which may be extracted by fitting a uniform stellar disk model to the visibility 
measurements. However, stars are not uniform disks, and one such deviation from a uniform 
disk is limb darkening, which is the apparent decrease in surface brightness towards the edge 
of the star. An observer will see deeper into a partially absorbing atmosphere when viewing 
the star at its center than when observing the limb, and deeper layers are hotter (brighter) 
than outer layers. This effect can be modeled, and the correction can be applied to diameter 
measurements for a more precise determination of the effective temperature (Code et al., 1976). 
A model-independent measurement of limb-darkening can be accomplished by analyzing visibility 
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data beyond the first lobe if visibility data have a high enough signal-to-noise. As a result, the 
uncertainty in the effective temperature is no longer limited by interferometric measurements, 
but rather by photometric measurements. 

Other applications related to angular diameter measurements arc the indirect calculation of 
stellar ages through isochrone fitting (Boyajian ct al., 2011). Asteroseismology and the study 
of stellar pulsation modes (Cunha et al., 2007) will soon become possible with long baseline 
interferometry due to the subpercent accuracy of angular diameter measurements. 

2.2 Fast rotating stars 

Fast rotating stars are particularly interesting targets for SII, since they are normally hot. 
Rapidly rotating stars are typically young stars of spectral types O, B, and A; some are indeed 
rotating so fast that the effective gravity in their equatorial regions becomes very small (at critical 
rotation even approaching zero), and easily enables mass loss or the formation of circumstellar 
disks. Rapid rotation causes the star itself to become oblate, and induces gravity darkening. 
A theorem by von Zeipel (1924) states that the radiative flux in a uniformly rotating star is 
proportional to the local effective gravity and implies that equatorial regions are dimmer than the 
poles. Spectral-line broadening reveals quite a number of early-type stars as rapid rotators and 
their surface distortion was already studied with the Narrabri interferometer, but not identified 
due to then insufficient signal-to-noise levels (Hanbury Brown et al., 1967). Clearly, high angular 
resolution images will enable to sec many of these interesting phenomena. 

A number of these fast rotators have now been studied with amplitude interferometers. By 
measuring diameters at different position angles, the rotationally flattened shapes of the stellar 
disks are determined. For some stars, also their asymmetric brightness distribution across the 
surface is seen, confirming the expected gravitational darkening and yielding the inclination of the 
rotational axes. Aperture synthesis has permitted the reconstruction of images using baselines 
up to some 300 m, corresponding to resolutions of 0.5 mas in the near-infrared H-band around 
A = 1.7 (Zhao et al., 2009a). In Figure 2.2, an image reconstruction of a-Cephei obtained by 
the Center for High Angular Resolution Astronomy (CHARA) is shown (Zhao et al., 2009b). In 
this image, we can clearly see the oblateness and pole brightening. Another interesting feature 
is that the bright spot does not appear to be exactly at the pole due to limb-darkening. 

With the relatively few stars that have been studied at high angular resolution, we have 
already been surprised at least once in the case of Vega. This star has always been one of the 
standard calibration stars for the apparent visual magnitude scale. When Vega was observed 
with the NPOI interferometer, it was learned that it is actually a fast rotating star viewed 
pole-on (Peterson et al., 2006). 
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Figure 2.2. Surface of a-Cephei. Stellar oblateness and pole brightenning can be seen. Courtesy 
of John Monnier. 

Predicted classes of not yet observed stars are those that are rotating both rapidly and 
differentially, i.e., with different angular velocities at different depths or latitudes. Such stars 
could take on weird shapes, midway between a donut and a sphere (MacGregor et al., 2009). 
There exist quite a number of hot rapid rotators with diameters of 1 mas or less. In fact, most 
hot (T > 10^° K) stars in the JMMC stellar diameter catalog have diameters smaller than 1 mas 
(402 out of 418 hot stars). Clearly the angular resolution required to reveal such stellar shapes 
would be 0.1 mas or better, requiring kilometric-scale interferometry for observations around 
A = 400 nm. 

A particularly interesting type of hot and fast rotating stars are Be stars. These are 
rotating at near their critical velocities and have strong emission lines and infrared excess, 
which provide evidence for the presence of a circumstellar envelope due to mass loss. The 
detection of partially polarized (~ 1%) light also indicates the presence of a circumstellar disk 
in most of these stars (Meilland et al., 2012). The study of the kinematics of the circumstellar 
material permit to further understand the nature of the mass loss of these objects. For example, if 
circumstellar matter is radiatively driven, then angular momentum conservation predicts that the 
tangential velocity scales as r^^. On the other hand, if the circumstellar disks are Keplerian, i.e., 
mechanically or viscosity driven, then tangential velocities scale as r^^/^. Such kinematic studies 
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Figure 2.3. a — Arae as imaged with the VLTI interferometer. 



can be done by noting that, for example, a fast rotating star viewed from its equator appears 
more red at one side, and more blue at the other. Therefore, images at different wavelengths 
reveal a shift in the image photo-center. Such studies are performed with the techniques of 
spectro-interferometry and spectro-astrometry (Oudmaijer et al., 2008), and have revealed the 
presence of both an equatorial disk which is mechanically driven (Oudmaijer et al., 2011), and 
a polar stellar wind, which is radiatively driven (Meilland et al., 2012). Such is the case of 
a — Arae, which was imaged with the VLTI array at 2.15 /xm (Meilland et al., 2007) as can be 
seen in Figure 2.3. 
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2.3 Radiatively driven stellar mass loss 

Another use of imaging in eonnection witli the astrophysics of hot stars is to quantify 
radiatively driven stehar mass loss. In radiatively driven stellar mass loss, matter is accelerated 
due to line scattering, which is the interaction of photons whose energy matches a particular 
energy level spacing in an atom (Castor et al., 1975). One current method of measuring total 
stellar mass loss is by analysis of P-Cygni^ spectral line profiles, whose signature is an asymmetry 
due to a blue shifted absorption. The blue shifted absorption in the spectral line is in turn due 
to the Doppler effect, i.e., as matter is accelerated by radiation, photons are red-shifted in the 
accelerated reference frame, and in an observers reference frame, only shorter wavelengths can 
meet the energy threshold for line scattering to occur. By analysis of P-Cygni spectral line 
profiles, only total mass loss rates have been measured so far^ (Puis et al., 2008). With high 
angular resolution imaging, it will become possible to map out the distribution of mass loss 
across the stellar surface as we shall now describe. 

An interesting characteristic of mass loss in these types of objects is that radiative transfer 
provides a connection between the luminosity map and mass loss map in the star. Much of the 
theory of radiatively driven stellar mass loss was developed by (Castor et al., 1975), and the 
most important result that can be derived in connection to high angular resolution imaging is 
that the luminosity map L{9x, 9y) is related to the local mass loss rate M{9x, 9y) by a power law 
of the form L{9x,9y) oc M{9x,9y)^"' . That is, an image of a star that is losing mass radiatively, 
provides a way to measure the mass loss rate at each point in the star. The exponent a can be 
shown to be 2/3 for hydrogen atmospheres, and has small deviations from this value when other 
elements are present (Puis et al., 2008). Even though the work by (Castor et al., 1975) is still 
relevant today, and is appealing due to its simplicity, it gradually fails as the star's mass loss 
departs from a steady state and when winds are optically thick. For this reason, models that 
allow for departures from local thermodynamic equilibrium (LTE) have been developed (Hillier 
Sz Lanz, 2001; Aufdenberg et al., 2002a) and used to study total mass loss rates in bright stars. 

It would be very interesting to simulate the capabilities of lACT arrays for high angular 
resolution imaging of mass loss in hot stars, and in particular, to quantify the capabilities of 
imaging variations of mass loss rates across the stellar surface. Prom observations of variable 
features within P-cygni lines and spectropolarimetric studies, there is increasing indirect evidence 
for the existence of regions in the stellar surface which exhibit higher mass loss rates at scales 



^P-Cygni is a Be star is the Cygnus constellation. 

■^Radio (cm) wavelength observations can also yield mass loss rates. The ionized winds of hot stars are free-free 
emitters and the flux scales as the density in the wind Abbott et al. (1980). 
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comparable to the size of the star and with time scales of days (Davies et al., 2005; Hamann 
et al., 2008). The origin of these features is unknown, but possibly related to magnetic activity 
or stellar pulsation. 

Hot O and B type stars are known to have strong radiatively driven winds, and they happen 
to be ideal targets for SII due to their high spectral density. As an estimate of the number 
of sources that can be imaged, according to the JMMC stellar diameter catalog, ~ 400 hot 
(T > 10000° K) stars can be imaged in less than 10 hours of observation time with a future lACT 
array such as CTA. Even before investigating high angular resolution imaging capabilities, by 
studying the expected fidelity of the second lobe in the visibility function, we will already gain 
insight into the the capabilities of lACT arrays to measure mass loss. There are currently no 
published results which analyze second lobe data for hot O or B type stars. 

In order to observe a localized region that is, for example, twice as bright as the rest of the star 
at blue wavelengths, it must experience a higher mass loss rate by a factor of 3 approximately. 
By extending existing non-LTE models such as the CoMoving Frame GENeral code (CMFGEN) 
Hillier Sz Lanz (2001), we should be able to predict brightness contrast values more precisely. In 
view of some preliminary results obtained here, stellar mass loss maps can very likely be imaged. 
The extent to which stellar mass loss can be imaged needs to be further investigated if we wish 
to test stellar atmosphere models. 

2.4 Binary systems 

From the list of fundamental stellar parameters that was mentioned in Section 2.1, a method 
for determining the stellar mass, arguably the most important stellar parameter, was not dis- 
cussed. The mass of a star essentially determines its fate. One way of finding stellar masses 
is through the determination of the orbital parameters of noninteracting binary systems. A 
considerable portion of the stars in our galaxy are in binary systems, so there is a large sample 
available to measure. In the case of noninteracting binary systems, the mass of the components 
does not change with time, so the determination of the mass allows us to test for stellar 
evolutionary models in general, i.e., not necessarily for stars belonging to multiple systems. 

When binary stars can be resolved, the orbital parameters can be found with a combination 
of spectroscopic and astrometric measurements. The spectroscopic data permit the evaluation 
of the orbital velocities along the line of sight, whereas the astrometric (and parallax) data allow 
us to evaluate the orbital velocity in the plane perpendicular to the line of sight. However, in 
order to accurately measure the orbital parameters, a considerable portion of the orbit has to 
be observed, and the number of observable binaries is reduced to those having orbits observable 
within human time-scales. Binaries with short orbital periods have small angular sizes and are 
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typically not resolved, so they have to be studied with long baseline interferometry techniques. 
Some of the orbital parameter measurements with long baseline techniques have been made 
with the Narrabri Stellar Intensity Interferometer (NSII), and more recently with the Sydney 
University Stellar Interferometry (SUSI), and the Very Large Telescope Interferometer (VLTI). 
So far, interferometric studies have measured masses in ~ 15 binaries with accuracies as small 
as a few percent (Quirrenbach, 2001). 

As binary components get closer to each other, a whole new set of phenomena start to occur. 
Numerous stars in close binaries undergo interactions involving mass flow, mass transfer and 
emission of highly energetic radiation: indeed many of the bright and variable X-ray sources 
in the sky belong to that category. However, to be a realistic target for interferometry, and 
intensity interferometry in particular, they must also be optically bright, which typically means 
B-star systems (Dravins, 2012). 

There has been a recent interest in studying massive (~ IOM0) binaries across the whole 
electromagnetic spectrum. The classical Be phenomenon, that was discussed in the previous 
section, has been associated to binary systems with hot massive stars. As mentioned above. Be 
stars possess both a polar wind and an equatorial decretion disk, and the degree to which each 
of these appear varies from one Be star to another. Binarity is thought to play a nonnegligible 
role, especially in the formation and/or truncation of the stellar disk (Millour et al., 2012). 

One well-studied interacting and eclipsing binary is /3 Lyrae (m^ = 3.5). The system is 
observed close to edge-on and consists of a B7-typc, Roche-lobe filling and mass-losing primary, 
and an early B-type mass-gaining secondary. This secondary appears to be embedded in a thick 
accretion disk with a bipolar jet seen in emission lines, causing a light-scattering halo above its 
poles. The donor star was initially more massive than the secondary, but has now shrunk to 
about 3Mq , while the accreting star has reached some 13Mq . The continuing mass transfer 
causes the 13-day period to increase by about 20 seconds each year (Harmanec, 2002). 

The first near-infrared optical image of the interacting binary system /3-Lyrae was recently 
obtained by the CHARA group (Zhao et al., 2008). With baseUnes up to 330m, the CHARA 
interferometer enabled the reconstruction of images at 2.2 /xm and 1.6 /xm which resolve both 
the donor star and the thick disk surrounding the mass gainer, located 0.9 mas away. A 
reconstruction obtained by the CHARA array is shown in Figure 2.4. The donor star appears 
elongated, thus demonstrating the photospheric tidal distortion due to Roche-lobe filling. 

Another type of related objects are X-ray binaries. These systems typically consist of a 
donor star and an accreting compact object, and can emit radiation as energetic as a few TeV. 
The high energy radiation most likely originates from high accretion rates or shocks from to the 
interaction between the stellar and pulsar winds. In Section 5.7, we will discuss the case of the 
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binary LS I + 61°303, which consists of a hot Be star and a compact object. This system may 
be too faint to be detected with interferometers, but it has been actively studied in the high 
energy (TeV) range by the VERITAS experiment (Acciari et al., 2009). 

Here we have already started to quantify the capabilities of imaging binary systems with 
lACT arrays (Nunez et al., 2012b), and some detailed results are given in Section 6.5.3. Imaging 
the effects mentioned above will further our understanding of stellar evolution, the formation of 
compact objects, type la supernova, and even the formation of planetary systems around young 
stellar objects (Karovska, 2006). 
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Figure 2.4. First resolved near optical image of the binary system /3-lyrae obtained by the 
CHARA array. The angular separation between components is 0.9 mas 



CHAPTER 3 



INTENSITY INTERFEROMETRY AND 



LIGHT FLUCTUATIONS 



3.1 Statistics of photo-electron detections 



Understanding the statistics of photo-electron detections is of central importance in inten- 
sity interferometry, since it ultimately determines the sensitivity of an intensity interferometry 
experiment. Here, we shall follow the discussion of chapter 9 in Mandcl & Wolf (1995). 

The intensity of light is in general a fluctuating random variable. However, we shall first 
study the case in which the intensity I{r, t) of the electromagnetic field is not a random variable, 
e.g. an ideal laser. The probability of detecting a photon in a small time St is proportional to 
the light intensity. Then the probability of detecting n photons in a time interval T is Poisson 
distributed, i.e., 



where the quantity in square brackets is the average number of detected photo-electrons, and rj 
is a constant that characterizes the detector. 

For a realistic light source (not an ideal laser), the intensity is actually a random variable 
whose distribution depends on the nature of the light source. Eq. 3.1 is true for a single element 
of the ensemble of the intensity. Therefore dp(n,t,T)/dT is not a Poisson distribution, but an 
average over equations of the form of 3.1, that is 




(3.1) 




(3.2) 



(3.3) 



(3.4) 



where 
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ft+T 

f^ = V Iir,t')dt', (3.5) 

and dV{/-J,)/dfj, is the probability distribution for ^u. Therefore, in order to find the probabil- 
ity distribution dp{n,t,T)/dT, we need to find dV{p,)/dfj, corresponding to the source under 
consideration. 

3.1.1 The statistics of a thermal source with a slow detector 

To find dV{n)/d^, we can first find dV{I)/dI, the probability distribution for the intensity 
(number of photons before going through the Poisson detector) . In the case of a thermal source, 
consisting of many uncorrelated oscillators, the electric field {E oc \/7) is Gauss distributed 
because it is the sum of many independent random variates (central limit theorem). Since the 
electric field is a complex quantity E = x + iy = Ae^'^ , the probability distribution of the real 
and imaginary parts is the product of both distributions, i.e., 

dxdy 2TTa'^ 

This can also be expressed as a function of the magnitude and the phase by noting that 
dxdy = AdAd(j) 

d'^PjA </>) ^ _A_ A2/2a^ 

dAdcp l^a-^ ■ ^ ' 

Here the probability distribution for the phase is dp{(j))/d(p = l/(27r), and the distribution 
for the amplitude is 

dA ~ a'''' ■ 

Now the distribution as a function of the intensity is found by noting that dI/dA = 2 A and 
that 2cr2 =< I >, so that 

which is an exponential distribution. 

We can calculate dji/dl to find dP{iJ,)/diJ, to finally calculate the probability distribution 3.4 
for detecting n photo-electrons. In general, this probability distribution is not Poissonian. 

We now consider the limiting case when /(r, t) undergoes large variations in the time T, or 
equivalently, when the electronic resolution time is much larger than the coherence time Tc- Here 
jj, has negligible fiuctuations and can be approximated by 
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,i = ri{I{r,t))T. (3.10) 

This in turn implies that dV{iJ,)/diJ, ^ 5{ijl — {I{r,t))T) (Dirac distribution), so that the 
distribution dp{n,t,T)/dT is Poissonian and we have reproduced eq. 3.1, corresponding to the 
case of no fluctuations. As soon as the distribution d'P(//)/d// has a finite width, we start to see 
deviations from Poisson statistics. 

3.1.2 The variance of the super- Poisson distribution 

We now calculate the variance of the probability distribution given by equation 3.4. Following 
(Mandel &: Wolf, 1995), The average number of photons in a time interval T is 



^ dp{n,t,T) 

(n) = — (3.11) 

n=l 

= r (3.12) 

'^^^d^^ (3.13) 
diJ, 

ifi) . (3.14) 



Following similar arguments. 



/ 2\ 2 dp{n,t,T) 

(n ) = — (3.15) 



n=l 



= </^' + /x>, (3.16) 

so that the variance is 



(An^) = (r? + {nf-2n{n)^ (3.17) 

= (^2^^) + (^)2-2(/x)2 (3.18) 

= (^2^ + (^) - (3.19) 

= (ra) + (A/x2). (3.20) 

The fluctuations in the detected number of photons reflect the fluctuations in the light 
intensity integrated over the resolution time. In the case of the thermal source and a "slow" 
detector (T S> Tc, where Tc is the coherence time), the probability distribution for the integrated 
light intensity ji can be found by using the central limit theorem. The variance of the intensity 
can be calculated from eq. 3.9 as 
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(Al2) = (3.21) 

= {if . (3.23) 

Within the electronic time resolution T, there are Tc/T possible intensity values, where Tc 
is the coherence time. Therefore, from the central limit theorem, as the number of different 
intensity values (T/tc) increases with electronic resolution time, the variance of the distribution 
of the integrated light intensity decreases as Tc/T. The variance of the integrated light intensity 

is then 

(A/x2) = (M)2(r,/T). (3.24) 
An important result that we have derived for the case of a slow detector is that 



(An^) = (n) + {ny 



(3.25) 



Here it is clear that the statistics are no longer purely Poissonian because the variance is no 
longer equal to the mean. These statistics are commonly known as "super-Poisson" statistics, 
and we can see from eq. 3.25, that deviations from Poisson fluctuations can only be seen with 
detectors whose resolution time is not too far away from the coherence time. A photo-detector 
measures a Poisson part, whose fluctuations are described by the first term in eq. 3.25, usually 
called shot noise, and a part related to intensity fluctuations, described by the second term in 
eq. 3.25, usually called wave noise. 

3.1.3 A Monte-Carlo simulation of photon-electron detections 

In Section 3.1.1, we studied the limiting case for a thermal source and a slow detector. 
We shall now consider the case in which deviations from Poisson statistics start to become 
visible. When the probability distribution for the integrated light intensity (eq. 3.25) is 
inserted in eq. 3.4, it is no longer straight-forward to derive an analytical expression for the 
probability of detecting n photons in a time T. However, a Monte-Carlo simulation can be made 
computationally by generating random events with this probability distribution. We simulate 
a source that produces on average 1.5 x 10^ photo-electrons per second approximately. For a 
realistic photo-multiplier tube, the resolution time is taken to be T = 10~^ s. The coherence time 
for thermal light is typically much shorter than the resolution time, but for illustrative purposes, 
the probability distribution is shown for a coherence time of Tc = 5 x 10"^*^ s in Figure 3.1, so 
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that deviations from Poisson statistics can be more easily seen. Here we took an observation 
time of 5 X 10~^ s, so that there are 500 simulated measurements. The standard deviation of 
the distribution shown in Figure 3.1 is approximately 5 photons per resolution time, which is 
one more photon per resolution time than if the distribution were purely Poissonian. A more 
detailed simulation of an intensity interferometer is described in Rou (2012) 



The wave noise can, in principle, be measured with a single detector. However, realistic 
sources have a smaller coherence time than the one considered in Section 3.1.3, and can be 
as small as 10~-^^s. Nevertheless, one can still see the effect by measuring the correlation 
between neighboring photo-detectors as illustrated in Figure 3.2. The Poisson fluctuations are 
not correlated between detectors, but the fluctuations due to intensity variations are correlated 
when both detectors are located within the same coherence volume (see section 1.3.1). This is 
because within the coherence volume, photons are indistinguishable and only symmetric states 
can occur between them, that is, they are correlated by the Bose-Einstein distribution. For the 
purpose of astronomical intensity interferometry, the correlation can be understood classically 
in the sense that both detectors are being driven by the same wave, and the wave has a definite 
frequency and phase within the coherence volume. Therefore, the wave noise is correlated 
between detectors in the same coherence volume. 

There is additional information contained in intensity correlations as we shall now see. Inten- 
sity interferometry allows us to measure correlations of intensities between pairs of telescopes, 
averaged over the signal bandwidth (denoted as (.•■)). If we express the instantaneous intensity 
as / = I + AI, where AI = I{t + dt) — I{t), the measurable quantity, denoted as 7^^^ is therefore 



3.2 Intensity correlations 
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^1/2) 



(3.26) 



ih) ih) 

((/i + A7i)(/2 + AJ2)) 
(Ji + All) {h + Ah) 
{1J2 + IiAh + AIJ2 + AI1AI2) 



(3.27) 



(3.28) 




= 1 + 



(AJ1A/2) 
ih) ih) ■ 



(3.29) 



In the previous expression we have used the fact that (A/) = 0. Eq. 3.29 expresses the fact 
that we measure correlations of intensity fluctuations. Now we make the connection between 
7^^) and the degree of coherence 7^^-' = 7. Assuming that the electric fields are Gaussian random 
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Figure 3.1. Probability distribution of the number of photons per resolution time. This 
corresponds to a source for which an average of 15 photons per resolution time can be detected. 
Here the electronic resolution time is T = 10"^ s, and the coherence time is Tc = 5 x lO"^'^ s. The 
pure Poisson distribution is compared with a Gaussian fit, and deviations from Poisson statistics 
start to become visible. 
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Figure 3.2. Schematic of the principle of an intensity interferometry experiment. If both light 
collectors are within the coherence volume of light, then the current (intensity) fluctuations are 
correlated. 
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variates and making use of the Gaussian moment theorem^, we can also rewrite eq. 3.26 as 



^ = {E,Et) iE,E*,) ^'-''^ 



\{EiEl 



|2 



= l + M'- (3.32) 
A comparison with eq. 3.29 yields the following important result. 



(3.33) 



Since 7 is the (complex) Fourier transform of the radiance distribution of the object in 
the sky, measuring intensity correlations enables the squared modulus of this quantity to be 
measured (Hanbury Brown, 1974). Therefore, the phase of the Fourier transform is lost during 
the measurement process, and it must be recovered for model- independent imaging (chapter 4). 

3.3 Signal- to- noise in intensity interferometry 

From the discussion of equation 3.25 we learned that the variance of photon counts, or 
equivalently the intensity fluctuations, contain two contributions: Poissonian shot noise, and 
wave noise, and it is the latter contribution that is correlated between detectors. When the 
degree of coherence is maximum (I7P = 1), the signal-to-noise ratio per electronic resolution 
time T is then the ratio of the wave noise and the shot noise, i.e., 

SNR = nj;. (3.34) 

Here n is the rate of detected photons within a certain optical bandwidth Ai/, and Tc is the 
coherence time. To further develop the previous expression, we can write the rate n/T as 

S = / PLd-. (3.35) 
T J dvdT ' ^ ' 

so that for a small bandwidth Az/ = I/tc 

n cPn' 1 

We should also note that when observing during a time Tq, the signal-to-noise increases as 
■sJTqJT, so the signal-to-noise becomes 



^Gaussian variates have the property that all higher order correlations can be expressed in terms of second 

order correlations. That is, 

{ztiZh ■ ■ ■ ZiNZjlZj2 ■ ■ ■ Zjn) = En! pairings (^'l^jl) • • • i^mZjN) 
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SNR=%^%Kfr2 



(3.37) 



where A/ is the electronic bandwidth and the factor of 2 is due to Nyquist's theorem. Now the 
rate of detected photons centered around a particular frequency can be written in terms of the 
rate of photons per area N (centered around a particular wavelength), the detector area A, and 
the quantum efficiency a as 



It is important to emphasize that is a property of the light source, therefore the signal-to- 
noise depends on the brightness of the object at the observed wavelength. It is also important 
to note that the signal-to-noise is independent of the optical bandwidth. The sensitivity cannot 
increase indefinitely by increasing A, since at some point the light collector will start to be 
large enough to resolve the light source and add uncorrelated intensities to the signal, therefore 
canceling the effect we want to measure. The integration time as well cannot be beneficially 
increased indefinitely since the finite point spread function of the optics results in integrating 
background light. Integrating background light places a limit on the minimum brightness the 
source can have, but does not pose a serious limitation for making precision measurements on 
a source that is much brighter than the night sky background. The electronic time resolution 
can in principle be as high as possible, but at some point the detected intensity fluctuations 
will be so fast that high precision optics are needed, therefore introducing additional technical 
difficulties associated with Michelson interferometry. A detailed discussion on the sensitivity of 
a modern stellar intensity interferometer is presented in Section 6.2. 




(3.38) 



and the signal-to-noise for an arbitrary value of |7p is now (Hanbury Brown, 1974) 



SNR = 7V^a|7|VroA//2. 



(3.39) 



CHAPTER 4 



PHASE RECOVERY 
4.1 Alternatives for imaging 

The imaging problem in intensity in interferometry is then reduced to finding the phase of 
the complex degree of coherence. There are several alternatives: The first is to live with the fact 
that the phase is not directly measured, and recover images using parametric models. In many 
cases, even when the phase is partially known, data are fit to a parametric model and relevant 
physical quantities are extracted from the model. 

The second option is to measure third order correlations between intensity fluctuations, 
similar to what is done in amplitude interferometry with the phase closure technique, and also 
similar to what is done in speckle interferometry. The phase problem in amplitude interferometry 
is discussed in section 4.3. 

The last option consists of tackling the phase retrieval problem from the Fourier magnitude 
data, and is discussed in sections 4.2 and 4.5. At first glance, the problem seems quite hopeless 
since any phase one postulates is consistent with the (measured) Fourier magnitude. However, 
we show in section 4.4, that since the Fourier transform of a function with finite support^ is 
analytic in the {u,v) plane, the phase can in principle be found by analytic continuation. We 
then present several approaches to the phase retrieval problem that make use of the theory of 
analytic functions in section 4.5. 

4.2 Phase retrieval problems in physics 

The phase retrieval problem and several related inverse problems arise in many fields of 
physics. Most of the phase retrieval problems arise when a wave is scattered off an object, then 
the information of the object is contained in both the magnitude and the phase of the propagating 
wave. When only the magnitude of the wave is measured, the phase is also needed to describe 
the object as accurately as possible. One of the earliest applications was in X-ray crystallography 
(Robertson, 1981), where a periodic crystal creates a difi"raction pattern corresponding to the 
squared modulus of the so called structure factor. The structure factor is the Fourier transform 
of the electron density function, and since only the squared modulus is measured, the phase needs 

^For example, a star a star 0{9) of angular size O has finite support, i.e., 0{6) = V 6 s.t. \9\ > O. 
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to be recovered to determine the crystal structure. Another example unrelated to astronomy and 
optics is found in quantum mechanical scattering (Russell et al., 1988), where the observable 
is the probability amplitude or squared-modulus of the wave function. The phase problem 
arises when one seeks the functional form of the interaction potential from knowledge of the 
squared-modulus of the wave function. This can be understood within the Born approximation, 
where the first order correction to the wave function of the outgoing particle is proportional to 
the Fourier transform of the potential. Only the magnitude of this Fourier transform can be 
measured, and the phase needs to be recovered to solve for the potential. Other examples where 
the phase problem is encountered include electron microscopy (Huiser &; Ferwerda, 1976) and 
wave-front sensing (Gonsalves &; Chidlaw, 1979). Here we concentrate on the application of the 
phase retrieval problem to astronomical intensity interferometry. 

4.3 Phase retrieval in amplitude interferometry 

To illustrate the fact that interference fringes contain phase information, consider the case 
of the following double slit experiment: Two narrow beams with a phase difference A between 
them interfere to form a diffraction pattern. Formally, the observed diffraction pattern is the 
intensity of the Fourier transform of the following radiance distribution 

B{e) = e'^/'^5{e + a/2) + e-'^l'^5{9 - a/2), (4.1) 
where the slit separation is a. The observed diffraction pattern is therefore 



= 2^|cos(fca;a + A)p, (4.3) 

where ^ is a constant specified by the detector characteristics. Therefore, the sinusoidal diffrac- 
tion pattern is displaced by and amount A. However, the problem is that besides the true phase 
difference A, there are also phase differences induced by atmospheric fluctuations in time-scales 
of the order of a few milliseconds (Labeyrie et al., 2006). Therefore, atmospheric fluctuations 
have the effect of drifting fringes in time. 

The approach used in amplitude interferometry is to apply a technique known as phase closure 
(Jennison, 1958). To briefly illustrate this approach, consider an array of detectors that can be 
divided into closed loops of three detectors (triangles) ijk. The signal at each detector contains 
an atmospheric phase shift (Ao,i, Aqj, Aq,*;). The phase at each detector is then cj)i = Ai — Ao,j, 
and the atmospherically modified coherence function between telescopes i and j is (Labeyrie 
et al., 2006) 
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lij=li3^^^''-^'-'\ (4.4) 

where = |7ij|e*(^'~^J^. Therefore, the product of the coherence functions along the closed 
loop is 

liil'iklki = l7i,-7infci|e^'+^^+^^ (4.5) 

The most important thing to note is that this quantity is independent of atmospheric turbu- 
lence as long as fringes are scanned in timescales smaller than atmospheric fluctuations (milli- 
seconds). The sum of phases in the exponential is a measurable quantity and is known as the 
closure phase. For an array of N receivers, there are N{N —l)/2 baselines and {N —1){N — 2)/2 
independent triangles. Therefore, in a nonredundant array, there are N — 1 more unknowns for 
the phase than there are closure phase equations. The procedure to find the phase consists in 
measuring the closure phase in each closed loop of the array,^ and the phase can be completely 
specified if there are enough redundant baselines. 

4.4 Uniqueness 

As was stated before, the phase of the Fourier transform has to be recovered from magnitude 
information only in intensity interferometry. To gain some intuition on the phase recovery 
problem, first consider the one-dimensional case of an object B{6) of finite extent. The Fourier 
transform of the one-dimensional object is an analytic function since it can be expressed as a 
z-transform, i.e., 

N 

j{z) = ^B{nAe)z''Ae, (4.6) 

n=0 

where z = exp (ik mAx AO) . 7(2;) is a polynomial in z and therefore an analytic function. An 
analytic function (of order zero) can in general be expressed as the product of its zeros, so eq. 
4.6 can be written as (Hadamard Factorization) 

N 

jiz) = cll{z-aj), (4.7) 

i 

so that all of the information of 7(2) is encoded in the roots aj. In SII we have knowledge of 
|7(z)p = 7(^)7(2""^), where 7(z) is a polynomial of degree N in z, and j{z~^) is a polynomial 
of degree N in z^^. The phase recovery problem can then be restated as finding 7(2;) from 
knowledge of 17(2:) p. The information contained in 7(2)7(2^"^) is also contained in 



^Two noncoUinear phase differences can be set to zero. 
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Q{z) = z''j{zMz-'), (4.8) 

which is a polynomial of degree 2N in z (Hayes theorem (Hurt, 1989)). The polynomial Q{z) 
has roots aj and aj^. The problem is then to find all the polynomials 7(2), with non-negative 
coeflttcients where either aj or aj^ is a root. If there are A'^ distinct roots, then there are no 
more than 2^ solutions. However, if there are N' roots that satisfy \aj\ = 1, then there are 
2N-N' golutions (Hurt, 1989)pg 30. For example, the Fourier transform of a step function 
(corresponding physically to a uniform disk-like star) has a corresponding z-transform of the 
form z^. The zeros of this function are all in the unit circle, so the solution is unique in this 
case. 

In general, given a polynomial 17(2;) p, the solution polynomial 7(2;) is not unique. However, 
all solutions are related to each other by a phase factor of the form^ Az^ {A, B e C). In image 
space, this is equivalent to solutions differing by translations and scale factors ((Klibanov et al., 
1995)). The set of solutions describing "the same object" are usually known as trivial associates. 

The statement of analyticity of the Fourier transform of a one-dimensional object is actually 
much more general than the discrete case treated so far. The Paley-Weiner theorem (Hurt, 
1989) states that the Fourier transform of a one-dimensional function with bounded support is 
an analytic function. The proofs of uniqueness are ultimately based on the uniqueness of analytic 
continuation. That is, if we have knowledge of a function in a region of the complex plane, by 
analytic continuation we can have knowledge of the function in the entire complex plane. ^ 

In the case of a two-dimensional function with compact support, its Fourier transform 
F{zx,Zy) is fully analytic (see Planchcrcl-Polya theorem). In two dimensions, an analytic 
function can in principle be factorized (Osgood product) in a similar form as eq. 4.7, but 
the form of each factor and the number of factors is unknown in general (Hurt, 1989). Zeros 
in two dimensions are not isolated, but rather form "lines" that uniquely define the function. 
This can be contrasted with the one-dimensional case, where the number of factors is known 
(eq. 4.7), and each factor corresponds to a root of the polynomial. The most important idea 
concerning uniqueness in two dimensions is that if F{zx,Zy) is irreducible, or cannot be written as 
the product of two analytic functions, then it is uniquely determined by \F{zx,Zy)\ (Sanz-Huan 
theorem). Going back to the discrete (polynomial) case, it has been noted (Hurt, 1989) that 
most two-dimensional polynomials are irreducible, and therefore uniquely determined up to 

^Provided there are no zeros in the origin. 

*A known example of analytic continuation is found in classical electrodynamics, when we wish to find the real 

part of the complex index of refraction, with knowledge of the imaginary part. These two quantities are related 
to each other by the Kramers-Kronig relations, also known as the Hilbert transforms. 
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trivial associates by their modulus. Moreover, even in cases when a polynomial F(zx,Zy) is 
irreducible, there is always a sufficiently small region around, say (-23;, 0; Zy,o) j where there exists 
an irreducible polynomial. 



The preceding section suggests that finding the zeros of the Fourier transform is a way of 
finding the solution. However, this a very unstable way of finding the solution (Klibanov et al., 
1995). Nevertheless, we shall see that care should be taken when extracting information of 
the phase when close to zeros of the Fourier transform. We shall now briefly describe some 
approaches to phase retrieval. 



Explicit formulae for the phase rely on the theory of analytic functions. The so called 
dispersion relations in particular are derived from the Cauchy integral formula by promoting the 
position variable x to be complex (Klibanov et al., 1995). Suppose that the degree of correlation 
as a function of position x can be expressed as 



For the moment, we assume that this function does not contain any zeros. Now we take the 
log of 7(x) and use the Cauchy integral formula^. The Cauchy integral along a large semicircle 
in the upper half complex plane is 



These integrals can be further simplified by making plausible assumptions of the asymptotic 
behavior of the magnitude and phase as the radius R of the semicircular path tends to infinity. 
For example, the magnitude can be assumed to decrease as l/x" for some n G M, and the phase 
can be assumed to be linear at infinity. The first term in the previous equation then becomes an 
integral along the real axis, and the second term becomes an angular integral on the semicircle 
that will not depend on x since R tends to infinity. Taking the real part of the previous equation 
yields 



4.5 Common approaches to phase retrieval 



4.5.1 Dispersion relations 



7(x) = |7(a;)|e^'^(^). 



(4.9) 




(4.10) 



5 
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P loQh(x)\ 

(j)(x) = — ; ^ ^^ dx' + ax + constant, (4.11) 

where P denotes the Cauchy principal value. The previous equation is also known as a logarith- 
mic Hilbert transform. When j{x) does contain zeros, there are additional contributions to the 
phase known as the "Blaske phase" 



X — a,- 
J ■' 



where aj refers to the zeros of j{x), and the problem is again reduced to finding the zeros of 
j{x). 

4.5.2 Cauchy-Riemann phase recovery 

Most of my phase retrieval research has concentrated in this method, which relies only on 
the theory of analytic functions, and which does not reduce to finding the zeroes of 7. We shall 
first study the one-dimensional case (Nunez et al., 2012b; Holmes Sz Belen'kii, 2004) and then 
treat the two-dimensional case to be used in SII analysis (Nunez et al., 2012b). 

4.5.3 The one-dimensional case 

If we denote I{z) = R{z)e^^^^\ where z = ^ + iip, we obtain the following relations from the 
Cauchy-Riemann equations^: 

dlnR ds , 
5$ _ aini^ _ ds 

'di ~ dip ~ W ^ ■ ^ 

where we have defined s as the log-magnitude. Notice the relation between the magnitude and 
the phase. By using the Cauchy-Riemann equations we can write the log-magnitude differences 
along the real and imaginary axes as: 



ds 5$ 

(4.17) 

If the log-magnitude were available along purely the ^ or the ip axes, we could solve the 
previous two equations for the phase. 



The C-R equations can be applied because "I" is a polynomial in z. 
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However, notice that because \z\ = 1, we can only measure the log-magnitude on the unit 
circle in the complex space 

In general, we can write the log-magnitude differences along the unit circle as 



Here A^± corresponds to phase differences along a direction perpendicular to As||, that is, 
perpendicular to the unit circle in the ^ — plane. We are however interested in obtaining A$||, 
so that we can integrate along the unit circle. 

The general form of $ can be found by taking second derivatives in eq. (4.14) and thus 
noting that $ is a solution of the Laplace equation in the complex plane. 

The general solution of ^{z) in polar coordinates {p,(p) is (Jackson, 1998) 

$(y9, (p) = ao + bo(f) + ^ (oj cos + bj sin j^) , (4.21) 
j 

where terms singular at the origin (p^^) have been omitted. These singular terms lead to 
ambiguous reconstructions including flipped images and have not been found to be essential for 
most reconstructions. 

Now taking the difference of $ along the radial direction we obtain 

A^^ip, (/)) = ^p^ ((l + - 1) (a,- cos + bj sin j<^) . (4.22) 

j ^ 

Note from eq. (4.19) that the length in the complex plane associated with As|| is A^ = 
|A^ + iAtpl, and that the length associated with A<I>_|_ is Ap = \ A^ + iAip\. Now setting p = 1, 
Ap = A(f), and for simplicity of presentation, expanding for small A^, we obtain 

j 

So now the coefficients aj and bj can be found using equations (4.18-4.19) from the measured 
As||, and thus $ can be found in the complex plane, with an uncertainty in oq and bo- The 
coefficients aj and bj can be calculated by performing the following integrals: 
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27rj 



cos j(f)d(t) (4.24) 



1 /■^'^ d$ 



bj = i^l ^smj(t>dct> (4.25) 
Note however that the previous expressions must exist, which is not the general case. More 
exphcitly, if the magnitude is zero, the log-magnitude is singular. When imaging finite objects 
in image space, there will always be zeros in the magnitude of the Fourier transform. In practice 
we are always sample limited and nothing prevents us from calculating aj and hj approximately. 

4.5.4 One-dimensional examples 

To illustrate the performance of the Cauchy-Ricmann phase reconstruction, some 
one-dimensional image reconstructions arc shown below. These examples do not include noise 
or realistic sampling of data. In Figure 4.1, the magnitude, phase, and reconstruction of a 
random image are shown. It should be emphasized that the only input in this example is the 
Fourier magnitude, and no prior information of the image for the reconstruction. A simpler 
example of a top-hat image reconstruction is shown in Figure 4.2. The main limitation of the 
Cauchy-Riemann algorithm in 1-dimension is due to the fit of eq. 4.23 by using 4.24 and 4.25, 
which results in not accurately reproducing discontinuities in the phase. 

4.5.5 Two-dimensional case 

Wc can think of this one-dimensional reconstruction as a phase estimation along a single 
slice in the Fourier plane. A generalization to two dimensions can be made by following the 
same procedure for several slices as described in Figure 4.3. In fact, the requirement that a 
two-dimensional complex function {zx,Zy) be analytic, is equivalent to satisfying the Cauchy- 
Riemann equations in both and Zy (Hormander, 1966). The direction of the slices is arbitrary, 
however for simplicity we reconstruct the phase along an arbitrary set of perpendicular directions 
in the Fourier plane, and noting that one can relate all slices through a single orthogonal slice, 
i.e., once the phase at the origin is set to zero, the single orthogonal slice sets the initial values 
for the rest of the slices. 

One can also require that the phase at a particular point in the complex plane be exactly 
equal when reconstructed along z^ or Zy since each reconstruction is arbitrary up to a constant 
(piston) and a linear term (tip/tilt). However, imposing this requirement results in a severely 
over-determined linear system. More precisely, by imposing equality in points in the complex 
plane, and having 2n slices (each with an unknown constant and linear term) , results in a linear 
system of equations and 4n unknowns. Alternative methods of requiring slice consistency are 
a possible way of improving phase reconstruction, but are beyond the scope of this work. 
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Figure 4.1. Example reconstruction of a random one-dimensional image. The top figure is the 
magnitude of the Fourier transform of the original image. The middle figure is the phase of the 
original image compared with the reconstructed phase using the Cauchy-Riemann algorithm. 
The bottom figure (in arbitrary units of intensity) is the original image and the image using the 
estimated phase. 




Figure 4.2. Example reconstruction of a top-hat function. The top figure displays the real 
and reconstructed phase using the Cauchy-Riemann phase reconstruction. The bottom figure 
displays the real and reconstructed image. 
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The Cauchy-Riemann approach, with horizontal or vertical slices, and a single orthogonal 

slice, gives reasonably good results; however, it is not the only possible approach. We have 
also investigated Gerchberg-Saxton phase retrieval. Generalized Expectation Maximization, and 
other variants of the Cauchy-Riemann approach. It is premature to conclude which of these 
approaches is best at this time, given the limited imagery and SNR levels that have been explored. 
However, the Cauchy-Riemann approach has shown to give better results in a number of cases 
(Holmes et al., 2010). 

4.5.6 Two-dimensional examples 

A few examples of two-dimensional image reconstructions are shown. Each of these examples 
takes the magnitude of the Fourier transform as the only input. In Figure 4.4, the reconstruction 
of an oblate object, e.g. a fast rotating star is, shown. In Figure 4.5, the reconstruction of 
a simulated image of the binary j5-lyrae is shown. As a final example, an image of Saturn 
is reconstructed in Figure 4.6. From the examples it can be seen that several main features 
are reconstructed approximately, and the quality of the reconstruction degrades with image 
complexity. More realistic examples are given in Chapter 6, as well as a more quantitative 
analysis of the reconstruction capabilities of this algorithm in the presence of noise, etc. 

4.5.7 Error-reduction algorithm 

The Gerchberg-Saxton algorithm, also known as the error-reduction algorithm, is an iterative 
procedure. Starting from a reasonable guess of the image whose phase is unknown, the algorithm 
consists in going back and forth between image and Fourier space, and each time imposing general 
restrictions. Since the data consist of Fourier magnitude measurements, the restriction in Fourier 
space is that the magnitude corresponds to the data. The restriction in image space can be as 
simple as requiring the image to be contained within some finite region. 

Figure 4.7 describes the Gerchberg-Saxton algorithm. Starting from an image Ofc, the first 
step consists in taking the Fourier transform to obtain something of the form M.^e^'^'-' . Now 
Fourier constraints can be applied, i.e., the magnitude is replaced by that given by the data, and 
the phase of the Fourier transform is kept. Now the inverse Fourier transform is applied and 
constraints can be imposed in image space. The constraints in image space can be very general, 
e.g. image positivity. However, if some apriori knowledge of the image is available, stronger 
constraints can be applied, and the algorithm converges faster. For example, if the image is 
known to have a finite size, a mask can be used, so that only pixels within the mask are allowed 
to have nonzero values. The performance of this algorithm depends strongly on the starting 
image, making it suitable for postprocessing purposes. Images using this algorithm are 
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Figure 4.3. Schematic representation of two-dimensional phase reconstruction approach. 
Several parallel slices are related to a single orthogonal slice. 
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Figure 4.4. Reconstruction of an oblate object, e.g., an oblate rotating star. Pristine image is 
shown in the bottom left corner 
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Figure 4.5. Reconstruction of a binary object. The pristine image (top left) is actually a 
simulated image of the well known binary system fi-lyrae. 




Figure 4.6. Reconstruction of Saturn. 
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Figure 4.7. Schematic of the Gerchberg-Saxton error reduction algorithm, 
presented in Chapter 7. 



4.6 Final remarks on phase recovery 

Phase retrieval is a field of research on its own right, and fully solving this mathematical 
problem has proven to be challenging indeed. The limitations of these algorithms, and reasons 
why some algorithms work better than others are still not fully understood (Hurt, 1989). 
However, it is clear that there is phase information contained in the Fourier magnitude, and 
perhaps one day, we will have full understanding of this mathematical problem. At this point, 
one is presented the following options: Either set on a quest to solve this problem, or use what 
is known so far to do science, e.g., astrophysics. I choose the later. The methods presented in 
the previous sections will be used Chapter 6 to quantify the imaging capabilities of future Air 
Cherenkov Telescope Arrays (lACT). 



CHAPTER 5 



AIR CHERENKOV TELESCOPE ARRAYS 
AND GAMMA-RAY ASTRONOMY 

Imaging Air Cherenkov Telescope arrays are primarily used for 7-ray astronomy, which 
investigates some of the most violent phenomena in the universe. In this chapter, the subject 
of 7-ray astronomy is briefly discussed. Even though the motivations for this field are entirely 
different from high angular resolution astronomy, they do share common interests for a few 
objects. One such object is the X-ray binary LSI + 61°303, which consists of a hot Be star, 
and a compact object, and was observed with the Very Energetic Radiation Imaging Telescope 
Array System (VERITAS). An analysis of 7-ray data allows us to constrain some fundamental 
parameters of the system (Nunez et al., 2011), and many remaining questions can potentially 
be answered with long baseline optical interferometry. 

5.1 Highest energy gamma-ray sources 

The earth's atmosphere is constantly being bombarded by very energetic charged particles 
known as cosmic rays, whose energy spectrum essentially follows a power law which spans 12 
orders of magnitude (10^ — 10^^ eV). Their origin is unknown since their angular distribution 
is isotropic, and questions such as acceleration mechanisms and energy dependent composition 
(e.g., single protons or heavy nuclei ) are still subject of debate. The field of 7-ray astronomy 
was initially proposed for finding the origin of cosmic rays. Photons are not deflected by the 
complex magnetic fields that isotropize cosmic ray detection, and studying the spectral energy 
distributions of photons helps determine the nature of the particle acceleration mechanisms. 

It has been 100 years since the discovery of cosmic rays, and their origin is still unknown, or 
at least highly debated. However, 7-ray astronomy is a flourishing field, and after the detection 
of the Crab Nebula as the first TeV 7-ray source in 1989, more than ~ 100 high energy (TeV) 
sources have been discovered. Figure 5.1 shows the sky map of 7-ray sources, which are divided 
in two main categories: galactic sources, which can be seen to lie along the galactic plane, and 
extra-galactic sources. Extra-galactic sources include active galactic nuclei (e.g., M87 Acciari 
et al. (2008a)), and more recently star-burst galaxies (e.g., M82 VERITAS Collaboration et al. 
(2009)). Galactic sources include supernova remnants, pulsar wind nebulae, X-ray binaries. 
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and unidentified objects. Rather than studying the possible 7— ray emission mechanisms, which 
include high accretion rates, or inverse Compton scattering from electrons accelerated by shock- 
waves generated by supernova explosions, we will analyze data from a particular high energy 
emitting binary system, and constrain some of its fundamental properties in section 5.7. 

5.2 Needs for gamma-ray astronomy 

The extremely small wavelengths associated to 7-rays {^^ 10~^^ m) do not allow for them to 
be detected with traditional optics such as mirrors since interactions are at the nuclear level. At 
these very high energies, large amounts of stopping material are needed, and this acts essentially 
as a calorimeter. In the case of GeV 7-rays, whose flux is of the order of the order of 10^^ cm s~^, 
enough material (~ Im'^) can fit in a satellite for them to be detected from space. Such is the 
case of the recent Fermi satellite, which has been extremely successful at detecting nearly 1000 
sources. 

As energies reach 1 TeV, the particle flux is of the order of 10~^^ cms~^, so very large areas 
(~ 100, 000 m^) are needed as well as vast amounts of stopping matter (1000 g/cm^), equivalent 
to 1 m of bricks! Detection from space becomes impractical, and in order to detect TeV 7-rays, 
the optically thick atmosphere is used to stop 7-rays, and large light (~ 100 m^) collectors detect 
the faint Cherenkov light produced as the electromagnetic particle showers propagate through 
the atmosphere. 

5.3 Imaging atmospheric Cherenkov technique 

When a 7-ray interacts with a nucleus at the top of the atmosphere, it induces an electro- 
magnetic cascade as illustrated in Figure 5.2. The interaction with the initial nucleus permits the 
creation of an electron-position pair, which then in turn emit radiation through Bremstrahlung 
when they encounter other charges. This process continues to develop and the shower continues 
to grow until particles reach an energy of a few hundred MeV and ionization dominates as an 
energy loss mechanism. At this point, e+e~ pairs are produced at a smaller rate, and the size 
of the electromagnetic shower starts to diminish. 

Since charged particles in the electromagnetic cascade travel faster than light in air, they 
emit Cherenkov radiation, analogous to the wake formed in water by a boat traveling faster 
than the speed of sound on the water surface. This light is seen as a "streak" of light in the 
focal plane of each telescope as shown in Figure 5.3. 

5.4 Analysis 

Most of the recorded data (99.9%) corresponds to cosmic ray induced showers, so much 
of the analysis has to do with discriminating 7-ray events from cosmic ray events. The main 




Figure 5.2. 



Schematic of electromagnetic cascade. 
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Figure 5.3. Image of shower in tiie camera of one telescope of the VERITAS array which 
consists of 4 telescopes. 
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difference between the two is due to the fundamental nature of the interaction (QED vs. QCD), 
and is reflected in the shape of the shower: cosmic-ray showers are typically "fatter" since pions 
generated in a collision have large transverse momenta compared with e+e^ pairs. In order to 
determine the initial energy of the photon, the analysis depends strongly on accurate shower 
simulations. By using the shower images in several telescopes, the source of 7-rays can be found 
through a geometric reconstruction. That is, a line is traced through the major axis of the 
"streak" seen in each camera. The intersection of these lines points to the source. 

Once 7-ray-like events are selected, the background needs to be subtracted. One way of 
accomplishing this is to first point the telescopes at the source, and then point away from the 
source, so that an estimate of the background can be found. This method is usually not practical 
since much time is spent looking at background. The way it is done in VERITAS is by pointing 
the telescopes at^ the source, and then selecting "off-regions" to estimate the background. 

5.5 X-ray binaries and 7-ray attenuation 

In the past few years, several high mass X-ray binaries have been detected as gamma 
ray emitters (Aharonian et al., 2006; Albert et al., 2009; Acciari et al., 2008b), causing an 
intensification of observational and theoretical interest. High energy emitting binary systems 
consisting of a main sequence star and a compact object are the only known variable galactic 
very high energy (VHE) sources, and their short periods of days or weeks make them even more 
interesting observational targets. These binary systems are starting to become astrophysical 
laboratories in the sense that by increasing spectral coverage and statistics, the nature of photon 
emission and absorption mechanisms is becoming more and more constrained. Here we will 
mainly be concerned with high energy (TeV) photons emitted from the vicinity of the compact 
object and interacting with the background black body radiation and ejected material from the 
companion star. Even though these systems can be incredibly complex, a simple model of the 
absorption mechanisms and how they affect the system's light curve can still shed light on many 
aspects such as the masses and the orbital parameters. 

One such example is the high energy emitting binary LS I + 61°303 (Massi et al., 2004). It 
was first detected in the TeV range with MAGIC (Albert et al. , 2009) and further observed with 
VERITAS at flux levels ranging between 5% and 20% of the Crab Nebula (Acciari et al., 2008b). 
This source has been observed throughout most of the electromagnetic spectrum starting with 
radio frequencies and extending to VHE gamma rays (Leahy, 2004). This broad spectral study 
indicates that the system consists of a main sequence Be star of mass Mi = 12.5 it 2.5 Mq 



^Telescopes are not actually pointed at the source, but observations are made in "wobble- mode" CITE. A 
detailed description, although interesting, is beyond the scope of this document. 
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(Casares et al., 2005), surrounded by a circumstellar disk (Grundstrom et al., 2007; Paredes 
et al., 2007), and a compact companion separated by tens of solar radii at periastron. The 
compact companion can be cither a neutron star or a black hole (Casares et al., 2005), and its 
exact nature is still subject of investigation and debate (Zdziarski et al., 2010). The maximum 
VHE emission occurs close to apastron (Acciari et al., 2008b; Albert et al., 2009), suggesting 
that absorption plays an important role in the modulation. 

In the following sections, we consider the attenuation of gamma rays due to interaction 
with background radiation and then consider the interaction with circumstellar material. The 
77 absorption mechanism in high energy emitting binaries was first pointed out by Gould &; 
Schreder (1967) and has been studied in the context of observed sources such as LS I + 61°303 
and LS'5039 (Dubus, 2006; Sierpowska-Bartosik Sz Torres, 2009). In the attenuation due to 
pair production, the two variables that play a main role are the concentration of background 
black body photons, and the energy threshold for pair production, which in turn depends on the 
scattering angle between the primary TeV photons and the low energy photons. 

5.6 Interaction with background radiation 

5.6.1 Radiative transfer equation 

In order to develop a gamma-ray attenuation model, we need to treat the general case 
of a binary system consisting of a VHE emitting compact object and a main sequence star. 
The radiative transfer equation (Chandrasekhar, 1960) for the intensity I{s,E), where s is the 
distance traveled by a photon of energy E from the emission point is 

^1^3. = _(i + cos n{s, e) a{E, e, E) + j{E, s) ; (5.1) 
as 

where n(s, e) is the spectral density of background photons of energy e emitted by the main 
sequence star, a{E,e,$,) is the cross section^ for the interaction between photons colliding at 
angle ^, and j{E, s) is a source term which we will now describe. 

5.6.2 Neglecting the source term 

Since the attenuation term is due to VHE photons creating e+e~ pairs, the source term is due 
to secondary gamma-rays in the electromagnetic cascade. The energy of these secondary gamma- 
rays is degraded by typically a factor of 4 with each interaction, and since VHE observations 
range between ~ 0.5 TeV to ~ 5 TeV (Acciari et al., 2008b), only those at the far end of the 
measured spectrum can contribute to the intensity at a fraction of their energy. However, as 



^Note that the term (1 + cos^(s')) simply corresponds to the relative velocity between the incident and target 

photons 



50 

we shall see in section 5.8.1, only photons in the lower part of the spectrum are attenuated 

considerably, and feed the development of the electromagnetic cascade. Photons in the far 
end of the observed TeV spectrum are considerably more scarce since the spectrum is steep. 
Consequently, we neglect the source term. 

Now we estimate the contribution of secondary gamma-rays with an over-simplistic model 
which helps justify our neglection of the source term in eq. 5.1. We can estimate the effect of 
secondaries as an increase in initial intensity /(sq + As,E) by 2/(so + As,4£^), i.e., instead of 
having I{so,E) in eq. 5.4 (defined below), we have 

/(so + As, E) /(so + As, E) + 2/(so + As, 4E)P{so + As, 4E), (5.2) 

where a photon of energy AE is assumed to produce an e'^e~ pair, which in turn emitt a 
gamma-ray of energy E. P(so -|- As, 4E) is the probability that the photon of energy 4E exists 
in the first place, and we have assumed an electromagnetic cascade toy model. Taking the 
intrinsic intensity to behave as a power law spectrum, /(so + As,E) = I{Eq) (^^^ ■, where 7 
is the spectral index, eq. 5.2 simplifies as 

/(so + As, E) /(so + As, E){1 + 2x 4r^P{sQ + As, 4E)). (5.3) 

Since 7 ~ 2 (see section 2.1) and P(so + As, AE) < 1, /(so + As, E) increases by a factor of 
~ 9/8 at most. 

5.6.3 Solution of the radiative transfer equation 

After neglecting the source term, the solution to the radiative transfer equation is 

/(s,E) = /(so,/^) expj-y ' {I + cos ^{s'))n{s',e')a{E,e',s')ds'de'^ . (5.4) 

Here sq is the emission point at the vicinity of the compact object (see Figure 5.4), and e 
corresponds to the threshold energy for pair production, 

_ mlc^ 

E{i+ cos i{s)r ^^-^^ 

Note that the dependence of the scattering angle ^ in eq. 5.4 has been changed to a depen- 
dence on the path s. The problem then reduces to calculating the integral in the exponential of 
eq. 5.4, also known as the optical depth t(s, E) (Rybicki & Lightman, 1979). In our calculation, 
we consider the main sequence star as point source, and in view of the results obtained by Dubus 
(2006), including the angular extension does not change our results significantly. 

The distribution of background black body photons can simply be taken as 
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n(r,S) = no(S)^, (5.6) 

where tq and no are the radius of the Be star and the density of background photons at this 
radius, i.e., 

Here, the photon density has already been integrated over the half sphere (solid angle). 

5.7 The case of LS 1+61 303 
5.7.1 Attenuation 

There is debate as to what is the mechanism responsible for high energy emission. However, 
the aim of this paper is not to model the gamma ray emission but rather to investigate the 
effects of attenuation. This allows us to derive a few characteristics of the main sequence star 
environment and compact object orbit. 

Grundstrom et al. (2007) reported a temperature of T ~ 3 x 10^ K and radius of i? ~ 6.7Rq 
for the Be star. The black body distribution peaks at a few eV, and the threshold energy for pair 
production with TeV incident photons is of the order of 1 eV, so that most of the background 
photons may contribute to the attenuation, provided the scattering angle is favorable. The 
background photon density is found to be of the order of n-y ^ 10^^ cm"^ at a the radius of 
the star. The circumstellar disk has been observed by Waters et al. (1988) and by Paredes 
et al. (2007), who estimate the disc ion density to be Ue ~ lO^^cm"^ at one stellar radius. The 
cross section for pair production is of the order of cr^^ O.lor at the threshold energy. The 
cross section for interaction with hydrogen has a constant value of a^n ~ 2 x 10~^(7t above a 
few hundred MeV (Heitler, 1954; Aharonian, 2004). With these cross sections, a first estimate 
suggests that both interactions may result in comparable degrees of attenuation. However, 
there is a strong angular dependence in the 77 interaction, the extreme case being when the 
both photons are emitted in the same direction, a configuration in which there is no VHE 
attenuation. Also the threshold energy is much higher when the incident and target photons 
are nearly parallel, so fewer background photons contribute to attenuation. Consequently, 77 
attenuation may not have strong modulation as a function of the orbital phase when compared 
with the modulation produced by interactions with the circumstellar material. 

5.7.2 Orbital parameters of LS I + 61°303 
The orbital parameters for LS I + 61°303 are still subject of research (Aragona et al., 2009; 
Grundstrom et al., 2007; Casares et al., 2005) and are sketched in Figure 5.4. The measurable 
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Figure 5.4. Sketch of the orbital parameters of LS I +61°303. The arrow points to the observer. 

quantities of interest are: The period P, the angle between the major axis of the ellipse and 
the line of sight ip, the projected semimajor axis (oisini), corresponding to the ellipse of the 
Be star^, the eccentricity e, the phase at periastron ^o, and the mass function f{mi,m2), which 
depends on the period and the radial velocity and relates the masses of both objects and the 
inclination angle i. The most recent orbital solution has been obtained by Aragona et al. (2009), 
where P = 26.4960 d, ^ = 40.5 ± 5.7°, aisini = 8.64 ± 0.52 i?0, e = 0.54 ± 0.03, (po = 0.275 
and f {mi, 1712) = 0.0124 it 0.0022 Mq. It is important to remember that the value of the angle 
i depends on the mass of the compact object, and our results may be used to constrain this 
quantity. Since the mass of the compact object is a function of the inclination angle, we will 
take this mass to be a free parameter of the model. 

5.8 The integrated flux of LSI + 61 303° 

Following observations of LS'/+61°303 from 09/2006 to 02/2008, the VERITAS cohaboration 
reported power law spectrum = <I>o ( '*') with a spectral index of 7 = 2.4:±0.2stat^0.2syst 
at energies above ~ 0.5 TeV, and between phases cp = 0.6 and <j) = 0.8 (Acciari et al., 2008b). 
Observations at lower energies made by Fermi between 08/2008 and 03/2009, indicate that the 
spectral index does not change significantly as a function of the orbital phase (Abdo et al., 2009). 
Therefore, we assume a constant intrinsic^ spectrum as a function of the phase at TeV energies. 



''The projected semi-major axis of the eUipse described by the compact object is typicaUy labeled as ai sini. 
*By intrinsic we mean non attenuated by pair production. 
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The integrated flux is then 



poo 

m = / 



I{E,ct))dE = Fo — IiE,cf>)dE, 



(5.8) 



dEdtdA 



where Eq depends on the detection threshold energy of the detector, and Fq is a normahzation 
factor that is taken as a free parameter. 



Figure 5.5 shows the attenuation as a function of the orbital phase for several different 
energies for the case of the compact object having the canonical neutron star mass {i ^ 64° 
or M I.5M0). In Figure 5.5 we essentially reproduce one of the results obtained by Dubus, 
except that the orbital parameters used are the newer set obtained by Aragona et al. (2009). 
When only interactions with the background black body photons are taken into account, and 
the orbital plane is closer to being seen edge-on, the optical depth approaches a minimum when 
the compact object is close to the main sequence star. This is especially the case for very 
high inclination angles, corresponding to the mass of the compact object being small, and close 
to the Chandrasekhar mass. This behavior can be understood from the angular dependence 
of the threshold energy in addition to the relative velocity of the incident and target photons 
approaching a minimum. Also, at high energies, the cross section for pair creation decreases as 
the inverse square of the center of mass energy, decreasing the optical depth even more. That 
is, even though the total density of background photons increases (as 1/r^) when the compact 
object approaches the Be star, a combination of the previously mentioned factors dominates as 
can seen in Figure 5.5. 

Figure 5.6 shows the normalized integrated flux assuming different inclination angles and 
corresponding compact object masses. The VERITAS data shown in Figure 5.6 (Weinstein, 
2008) were binned to show a single light-curve as opposed to monthly data. If we assume that 
the emission comes from the vicinity of the compact object, and is isotropic, and constant as 
a function of the orbital phase, then these results lead us to conclude that there must be an 
additional attenuation mechanism at play. 

5.8.2 Light curve including 77 and interactions 

The detailed structure of the circumstellar material surrounding a Be star in the presence 
of the compact companion has been studied in detail by Waters et al. (1988), Marti &; Paredes 

(1995) 



5.8.1 Light curve assuming only 77 interactions 
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Figure 5.5. Attenuation e "^^^ as a function of the orbital piiase for different incident photon 
energies (77 interactions only). A mass of I.SMq was assumed for the compact object. 
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Figure 5.7. Light curve including an isotropic distribution of circumstellar material composed 
oh hydrogen. The mass of the compact object was set to A.QMq for the two curves whose peak is 
at phase ~ 0.7, and 1.5Mq for the curve whose peak is closer to phase ~ 0.4. The concentration 
hydrogen at r 2± lOOi?0 for each curve is labeled in the top right-hand corner of the figure. 



and Reig et al. (2000) among others. It is thought to have a main equatorial disk-like component 
and a polar wind. Typically, the parameters that describe the dccrction disk include: The mass 
loss rate, the wind termination velocity, the half opening angle of the disk, and the radius of 
the disk. When comparing the quality of the data shown in Figure 5.6 and the complexity 
of the models that describe the circumstellar material, only an order of magnitude estimate 
of the density of the material and its extension in the system can be achieved. With this in 
mind, we rather assume a simple isotropic distribution of material that decreases as a power q 
of the distance from the Be star (n = no(ro/r)*). We start by setting q = 2 and then consider 
different radial dependences for comparison. Parameters found from existing models are taken 
into consideration for our approximation. We assume that most of the material surrounding the 
Be star is composed of hydrogen, whose cross section with high energy photons is approximately 
(Heitler, 1954) a-yH — 2 x 10~^crr and roughly independent of the energy above a few hundred 
MeV. 

We can now add this contribution to the optical depth and obtain the light curves shown in 
Figure 5.7, and thus constrain the mass of the compact object and the density of the circumstellar 
material at a distance of r lOOi?© (characteristic order of magnitude of the system). For 
the case of a constant cross section and a distribution of hydrogen, the optical depth can 
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actually be found analytically as can be seen in the supplementary Section 5.12, and the behavior 
at the phase near emission peak is roughly Gaussian. To find the best fit values we take the mass, 
characteristic density and normalization factor as free parameters. As mentioned in the previous 
section, higher masses for the compact object shift the emission peak to higher values of the 
phase as shown in Figure 5.7. From this figure it is clear that the emission peak corresponding 
to a canonical 1.5Mq neutron star is only marginally supported by observations. 

We take the inclination angle, characteristic density, and normalization factor as free param- 
eters. We find the inclination angle to be i < 28° (M2 > 3Mq) at the 89% confidence level 
(CL) in the context of this model, and i < 34° (M2 > 2.5M0) at the 99% CL. These limits 
are not in good agreement with the neutron star scenario generally favored for the broad-band 
spectrum it implies^. However, our results are still consistent with other observational con- 
straints (10° < i < 60°) (Casares et al., 2005) obtained from optical spectroscopy. As for the 
circumstellar material, if we assume the characteristic extension to be ro ~ lOOi?0, consistent 
with more sophisticated models (Sierpowska-Bartosik Sz Torres, 2009), then the density of 
hydrogen in the disk is found to be 2.0 x lO^^cm"^ < uh < 1.9 x lO^^cm"^ at the 99% CL 
and UH = i^-^^Y) X lO^^^cm-^ at the 68% CL. 

By integrating the volume density along the line of sight to the compact object at apastron, 
we find a column density of 1.9 x lO^^cm-^ < Nr < 1.8 x lO^^cm-^ at the 99% CL, which 
is much higher than results found elsewhere in the literature (Waters et al., 1988; Marti Sz 
Paredes, 1995; Esposito et al., 2007). In particular, when we use the column density found by 
X-ray observations Nh = (5.7 it 0.3) x lO^^cm"^ Esposito et al. (2007), we find a reduced of 
3.06 (11 degrees of freedom), corresponding to a probability Pix^ > 3.06) = 0.04%. A rough 
estimate suggests that by including 10% of helium, the column density would be reduced by 
a factor of ~ 2, which is not sufficient to achieve compatibility with X-ray results. 

Density profiles in Be stars typically have radial dependences of 1/r'^, where 2.3 < q < 3.3 
(Lamers & Waters, 1987), depending on the opening angle of the disk. Therefore we expect our 
constraint on the density to constitute a lower bound^. We perform our calculation with q = 3 
and note that our results do not change considerably. 

The hydrogen density also corresponds to a mass loss rate of Mi ^ 1Q~^U Meyr"^, 
where Q, is the solid angle. Typically accepted values for the mass loss rate are in the range of 
~ lO^^Mgyr-^ to lO^^Moyr^^ as have been reported by Snow (1981) and Waters et al. (1988) 
among others. A first glance at our result for the mass loss implies that it does not agree with 

^See Zdziarski et al. (2010) for more details 

®This is assuming tliat the disk and orbit lie in the same plane. 
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the observations, i.e., setting = 47r and Vwind ^ 100km s^"*^ (Waters et al., 1988). However, if 
we relax the assumption of an isotropic distribution of hydrogen, our result implies that small 
solid angles are favored as well as small velocities for the stellar wind. Small solid angles are 
consistent with the thin disk scenario that is most commonly accepted. Small velocities of the 
order of a few kms~^ are however not consistent with what is found elsewhere in the literature, 
e.g., Waters et al. (1988), and the wind indeed has higher velocities, this would imply that the 
system may have been observed while in a state of high mass loss rate. 

5.9 Discussion on LSI + 61 303° 

Since, in the TeV range, the interaction with matter is approximately independent of the 
energy, and since, as Figure 5.6 shows, 77 interactions are insufficient to account for the orbital 
modulation, then the intrinsic nonattenuated differential spectrum is essentially the same as the 
observed spectrum (a power law of spectral index —2.4). However, the intrinsic TeV luminosity 
is several orders of magnitude higher than the measured luminosity. Taking the distance to the 
source to be approximately 1.8 kpc (Steele et al., 1998), we find the intrinsic luminosity to be 
L 5 X lO^^ergs"^ when the hydrogen density is of the order of ~ 5 x lO^^cm"^. This intrinsic 
luminosity is comparable to that suggested by Bottcher (2007) for LS 5039, the only other known 
TeV binary thought to contain a black hole. 

It is interesting to compare this intrinsic luminosity to the Eddington luminosity^ LEdd ~ 
1.3 X 1O^^(M2/M0)ergs~"'^, which is comparable to L, and implies that radiation may be beamed 
in our direction. It is also interesting to calculate the accretion rate that would be needed in 
order to obtain the intrinsic luminosity: By taking L ~ GM2M2/-R, where R is of the order of the 
Schwarzschild radius (2GM2 /(?)■, we find M2 ~ 2 x lO^^MQyr"-^. This rate is comparable with 
the observed mass loss rate of ^ lO~^M0yr~^. The fact that the accretion rate is comparable 
to the measured mass loss rate, suggests that the fiow of matter can be quite complicated, e.g., 
an increase in the accretion rate would strip most of the circumstellar mass, leading to time 
variability. This may explain the fact that no VHE detections have been reported since 2008. 

Still assuming the intrinsic luminosity to be constant in time, we can estimate the amount 
of hydrogen needed to attenuate the source to below the detectability threshold. We find that 
the density must increase from ~ 5 x lO^^cm"^ to ~ 5 x lO^^cm"^ at the characteristic distance 
of lOOi?0. This amount of hydrogen in turn leads to much higher mass loss rates than those 
observed, and it may also imply a stronger activity of the source. 

It is worth mentioning that the attenuation model is not the only possible way to account 



At the energies considered here, the cross section for inverse Compton is ~ O.Iut 
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for the modulation. For example, there is also the possibility of the emission being anisotropic, 
and the modulation resulting from a geometrical effect. This possibility is described in detail by 
Zdziarski et al. (2010), were a shocked pulsar wind with a large Lorentz factor is thought to be 
the cause of emission. 

5.10 Final remarks on LSI + 61 303° 's 
TeV data analysis 

For the case of LS I + 61°303, we find that attenuation due to 77 interactions with the 
background radiation docs not account for the observed high energy flux modulation as a function 
of the orbital phase, namely a narrow peak near apastron. This effect leads us to investigate 
some properties of the ejected material from the Be star, and the inclination angle of the orbit. 
We find the angle of the orbit to be i < 34° (M > 2.5Mq) at the 99% confidence level, suggesting 
that the compact object is a black hole rather than a neutron star. We also find the density 
of hydrogen in the disk to be 2 x lO^^cm-^ < ur < 2 x lO^^cm'^ at the 99% CL (at WORq), 
which accounts for most of the observed gamma ray absorption. If the compact object is indeed 
a black hole as our analysis suggests, then the gamma ray emission is likely to be powered by 
accretion (Zdziarski et al., 2010). Also, a black hole scenario might be even more complicated 
due to the possibility of VHE emission originating from termination of jets, therefore we cannot 
exclude the possibility of the modulation being due to geometrical effects. Current VHE data 
do not allow us to constrain the system much more than what we have already done, and the 
fact that VHE detections have not been reported since the VERITAS (Acciari et al., 2008b) and 
MAGIC (Albert et al., 2009) detections where made, makes the problem even more puzzling. A 
possible explanation might originate from a complex matter flow. This is suggested by the fact 
that the accretion rate needed to explain an intrinsic nonattenuated luminosity, is comparable 
to the measured mass loss rate of the Be star. 

An inconsistency arises when comparing our results with those derived from X-ray observa- 
tions. We find the column density to be 1.9 x lO^'^cm'^ < Nr < 1.8 x lO^^cm"^ (-99% qL), which 
is only compatible with X-ray results at the 0.04% confidence level. Such an incompatibility may 
imply that pair production in the stellar wind is not the cause of the modulation. Consequently, 
our estimates on the mass and column density may not be valid. An alternative explanation 
by Zdziarski et al. (2010) suggests that the modulation is due to a geometrical effect. Here 
a shocked pulsar wind is thought to flow along a cone with a large Lorentz factor, producing 
beamed radiation which can be seen when the cone passes through the line of sight. 
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5.11 Prospects for LSI + 61 303° 
at high angular resolution 

Imaging at high angular resolution will allow us to further understand the nature of this 
object, and more objects of this class. The angular size of the Be star is approximately 
0.3 mas, and the angular separation between components may be larger by an order of mag- 
nitude depending on the inclination of the orbital plane and the orbital phase. Therefore, only 
interferometric techniques allow us to resolve this system. Radio observations with the Very 
Large Baseline Interferometer (VLBI) show structure at the milliarcsccond scale (Massi et al., 
2004) and show evidence of a precessing jet associated to the compact object. However, more 
information about the circumstellar environment of the Be star can be obtained by going to 
shorter near-infrared wavelengths since Be stars are known to have expanding dust shells, viscous 
disks, and/or strong radiatively driven winds. Current instruments such as CHARA, whose 
angular resolution can be as good as 0.3 mas at 550 nm, could use their largest 330 m baseline to 
obtain spectro-intereferometric data, where a shift in the image photocenter as a function of the 
wavelength may allow us to constrain the kinematics of the circumstellar matter^. If observations 
are done in the K band (~ 2200 nm), the angular resolution will decrease to ~ 1.3 mas, but it 
would be interesting to measure the interferometric visibility across the Hq, emission line, which 
is associated to the cool circumstellar environment. If a decrease in the visibility is evident, this 
would imply that the 330 mas baseline resolves the circumstellar environment, if no decrease is 
evident, an upper limit to the physical extension can be found. 

To obtain a fully reconstructed optical image, much better baseline {(u,v)) coverage is 
necessary, and going to shorter wavelengths may be beneficial in terms of angular resolution. At 
these short wavelengths, information of the stellar shape and temperature distribution can be 
obtained. In order to image features ranging between 0.3 — 3 mas, an instrument would require 
baselines ranging between a few tens of meters to a few hundred meters. In terms of angular 
resolution, this is within the capabilities of future intensity interferometers, whose simulated 
results show that imaging stellar shapes and temperature distributions is indeed possible (see 
Chapter 6). However, this object is just barely within the brightness detectability limit with 
intensity interferometry, and more detailed simulations are needed in order to determine if this 
is a suitable target. 



the Be star is observed edge-on, then one side should be blue shifted, and the other should be red-shifted 
since it is fast rotating. The measured phase of the complex visibility would be consistent with a nonsymmetric 
object. 
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5.12 Supplement: Optical depth for constant 
cross section and density distribution 

Using a distribution of hydrogen, the cross section accounting for interactions 

between VHE photons and hydrogen, and the system of coordinates shown in Figure 5.8 (corre- 
sponding to an orbital plane seen edge on) , we can calculate the integral for the optical depth 
to be 



f 

•JXi 



+ yf + zf 



dx 




(5.9) 



(5.10) 



where Xi, yi, and Zi are functions of the orbital angle 0, and tq is the characteristic radius of the 
hydrogen disk. For the case of a circular orbit as seen edge on (Figure 5.8), we can easily see the 
limiting behavior of the intensity as a function of the orbital angle. That is, expanding around 
^ ~ reveals that the attenuation around this region behaves like a Gaussian. 



ForOr^O: I{ri,e) = Io{e,ri)e 



(1+^2 



(5.11) 



Similarly, expanding around ~ 7r/2 reveals that the attenuation behaves like a decreasing 
exponential 



For 0^Tr/2: I{ri,9) = h{9,n)e 



(5.12) 



For a more complicated geometry of LS I + 61°303, it is now just a matter of inserting the 
appropriate expressions for Xi{9), yi{6) and Zi{9). 
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Figure 5.8. Upper figure: Coordinate system used for calculating the optical depth (eq. 5.10). 
Lower figure: Light curve for a circular orbit. 



CHAPTER 6 



AIR CHERENKOV TELESCOPE ARRAYS 
AS SII RECEIVERS 

6.1 A revival of SII 

Even though Stellar Intensity Interferometry (SII) was abandoned in the 1970s, there has 
been a recent interest in reviving this technique, mainly due to the unprecedented {u, v) plane 
coverage that future imaging air Cherenkov telescope (lACT) arrays will provide (Consortium, 
2010). The possibility of probing stars at the submilliarcsecond scale and visible wavelengths 
has motivated new developments in instrumentation and simulations, the latter being the focus 
of this chapter. 

Recent results obtained with amplitude (Michelson) interferometry have started to reveal 
stars as extended objects (e.g., Baldwin et al. 1996; Pedretti et al. 2009), and with nonuni- 
form light intensity distributions in the milliarcsecond scale. Such interesting results can be 
further investigated with SII taking advantage of the longer (km) baselines and relative ease of 
observing at shorter (blue) wavelengths. For example, measuring stellar diameters at different 
wavelengths, will make it possible to further investigate the wavelength dependence of limb 
darkening, (Mozurkewich ct al., 2003) and thus constrain stellar atmosphere models. Radii 
measurements with uncertainties of a few percent, along with spectroscopic measurements are 
necessary to constrain the position of stars in the HR diagram (e.g., Aufdenbcrg ct al. 2005). 
With the methods described in this chapter, we show that diameters can in principle be measured 
with accuracies better than 1% when using realistic array configurations for future experiments 
such as CTA (Cherenkov Telescope Array). As another example we can consider fast rotating 
B stars, which are ideal candidates for imaging oblateness, pole brightening (Monnier et al., 
2007; von Zeipel, 1924), radiatively driven mass loss (Friend & Abbott, 1986), and perhaps even 
pulsation modes (Saio et al., 2006). The impact of rotation on stellar evolution is nontrivial, 
and several studies have been made in the subject (e.g., Martin Sz Claret 1996; Maeder 1997). 
Images of rotating stars have become available in the past few years (e.g., Monnier et al. 2007; 
Aufdenberg et al. 2006), and measurements of oblateness with accuracies of a few percent have 
been made. We will show that this is comparable to what can be achieved with SII using large 
arrays of Cherenkov telescopes. There is also the case of interacting binaries, for which we 
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can measure angular separation, diameters, and relative brightness. It may even be possible 
to measure mass transfer (Verhoelst et al., 2007). Measurements of the angular separation in 
binaries is crucial for determining the masses of stars. These masses must be found to within 
~ 2% (Andersen, 1991) in order to test main sequence models. With the methods described 
in this chapter, wc show that angular separations can be found to within a few percent from 
reconstructed images. 

In preparation for a large-scale SII observatory deployment, several laboratory experiments 
are in progress (LeBohec et al., 2010). Their main goal is to measure light intensity correlation 
between two receivers. It is also worth mentioning the StarBase (LeBohec, 2007) observatory 
(located in Grantsville, Utah) which consists of two 3m light receivers separated by 24 m and 
which will be used to test high time resolution digital correlators, band to measure the second 
order degree of coherence for a few stars (see chapter 8). Various analog and digital correlator 
technologies (Dravins et al., 1994) are being implemented, and cross correlation of streams of 
photons with nanosecond-scale resolution has already been achieved. 

Intensity interferometry, unlike amplitude interferometry, relies on the correlation between 
intensity fluctuations averaged over the spectral band at electronic (nanosecond) time resolution. 
These averaged fluctuations are much slower than the (femtosecond) light wave period. This 
correlation is directly related to the complex degree of coherence 7^ as (Labeyrie et al., 2006) 

<h><Ij> ^ ' 

Here, < Ij > is the time average of the intensity received at a particular telescope i, and 

A/j is the intensity fluctuation. Measuring a second-order effect results in lower signal-to-noise 
ratio when compared to amplitude interferometry (Le Bohec & Holder, 2006). This sensitivity 
issue can be dealt with by using large light collection areas (such as those available with air 
Cherenkov telescopes), longer exposure times and baseline redundancy. 

The low frequency fluctuation can be interpreted classically as the beat formed by neighboring 
Fourier components. Since SII relies on low frequency fluctuations, which are typically several 
orders of magnitude smaller than the frequency of optical light, it does not rely on the phase 
difference between light waves, but rather in the difference between the relative phases of the 
two components at the detectors (Hanbury Brown, 1974). The main advantage is the relative 
insensitivity to atmospheric turbulence and the absence of a requirement for sub-wavelength 
precision in the optics and delay lines (Hanbury Brown, 1974). 

The complex mutual degree of coherence 7 is proportional to the Fourier transform of the 
radiance distribution of the object in the sky (Van Cittert-Zernike theorem). However, since 
with SII, the squared-modulus of 7 is the measurable quantity, the main disadvantage is that the 
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phase of the Fourier transform is lost in the measurement process. The loss of phase information 
poses a severe difficulty, and images have in the past been reconstructed from the bispectrum 
technique, using monolithic apertures (e.g., Lawrence et al. 1990). The imaging limitations can 
be overcome using a model-independent phase recovery technique. Even though several phase 
reconstruction techniques exist (Ficnup, 1981), we concentrate on a two dimensional version of 
the one dimensional analysis introduced by Holmes and Belen'kii (2004), which is based on the 
theory of analytic functions, and in particular the Cauchy-Riemann equations. 

Following recent successes in Gamma ray astronomy, a next generation Chcrenkov telescope 
array is in a preparatory stage. This project is currently known as CTA (Cherenkov telescope 
array) (Consortium, 2010), and will contain between 50 and 100 telescopes with apertures ranging 
between 5 m and 25 m. In this chapter we investigate the sensitivity and imaging capabilities of 
SII implemented on such an atmospheric Cherenkov telescope array. We start with a discussion of 
sensitivity (section 6.2), followed by a discussion of simulating noisy data as would be realistically 
obtained with such an array (section 6.3). Since data have a finite sampling in the {u,v) plane, 
we discuss our method of fitting an analytic function to the data in order to estimate derivatives 
which are needed for phase reconstruction (section 6.4.1). We then proceed to quantify the 
reconstruction quality using several criteria. We start with the simple case of uniform disks 
(section 6.5.1) and progressively increase the degree of image complexity by including oblateness 
(section 6.5.2), binary stars (section 6.5.3), and obscuring disks and spots (section 7.1). 

6.2 Sensitivity 

The signal to noise ratio (SNR) for an intensity correlation measurement depends on the 

degree of correlation 7, the area A of each of the light receivers, the spectral density n (number 
of photons per unit area per unit time, per frequency), the quantum efficiency a, the electronic 
bandwidth A/, and the observation time t. The SNR can be expressed as (Hanbury Brown, 
1974) 

SNR = n{X,T,my)AaJ^^/Aft/2. (6.2) 

This expression can be found from the results presented in section 3.1.2, and eq. 3.25 in 
particular. The SNR is essentially the ratio between the wave noise, which is correlated between 
neighboring detectors, and the shot noise, which is uncorrected between detectors. 

The spectral density n is related to the visual magnitude niy of the star as well as its 
temperature T and observing wavelength A. The spectral density n(A, T, rriy) is the number of 
black body photons per unit area, per unit frequency and per unit time. The dependence of 
the visual magnitude rriy is found by recalling that the flux for a O*'^ magnitude star with a 
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temperature of 9550°K observed at 550 nm is 3.64 x lO^^^Wm^^Hz^^ (Bessell, 1979). This in 
turn corresponds to a spectral density of 10^"^ m~^s~^Hz~^. The spectral density as a function of 
temperature (for different visual magnitudes and observing wavelengths) is shown in Figure 6.1, 
and we see that at constant visual magnitude and observing wavelength, higher temperatures 
correspond to higher spectral densities. We find that the increase in temperature AT(A, T, Amy) 
is approximately Am^ for the range of temperatures and wavelengths considered in Figure 
6.1. For example, at 400 nm, a decrease of 1 visual magnitude is equivalent to increasing the 
temperature of the star from 5000 K to 5700 K. Therefore, bright and hot targets are the most 
easily detected with SII. 

We use a preliminary design of the CTA project as an array configuration (Bernlohr, 2008), 
which is shown in Figure 6.2. This array contains iV = 97 telescopes and N{N — l)/2 = 4646 
baselines (many of which are redundant) which are shown in Figure 6.3. Each detector is 
assumed to have a light collecting area of 100 m^ and a light detection quantum efficiency of 
a = 0.3. Using a \/D criterion, we find that the largest baselines of 1.5 km resolve angular scales 
of ~ 0.05 mas at 400 nm. The smallest 48 baselines of 35 m resolve angular scales of ~ 2 mas . 
However, we show in section 6.5.1, that the largest angular scales that can be realistically imaged 
with our analysis, in a model independent way, are more determined by baselines of ^ 70 m. 
This is because the estimation of derivatives of the phase (needed for phase recovery) degrades 
as the number of baselines is reduced. Baselines of 70 m resolve angular scales of ~ 1.2 mas at 
400 nm. 

These order-of-magnitudc considerations are taken into account when performing simulations 
and image reconstructions, i.e., the minimum and maximum size of pristine images that can be 
reconstructed by data analysis, do not go far beyond these limits. More precise array resolution 
limits are presented in section 6.5.1 (diameters ranging between 0.06 mas — 1.2 mas) . By 
combining these angular scales with the SNR (eq. 6.2), we obtain Figure 6.4. This figure shows 
the highest visual magnitude, for which photon correlations (with I7I = 0.5) can be detected 
(5 standard deviations), as a function of the temperature, and for several different exposure 
times. Also shown in Figure 6.4, is the shaded region corresponding to angular diameters 
between 0.03 mas and 0.6 mas ^, and observable within 100 hrs. From the Figure we can see 
how correlations of photons from faint stars can be more easily detected if they are hot. To 
quantify the number of stars for which photon correlations can be detected with the lACT 



^ These curves of constant angular size can be found approximately by recalling that the visual magnitude m„ 
is related to a calibrator star of visual magnitude mo by: (m„ — mo) = — 2.51ogF/Fo. Here, F and Fg refer to 
the flux in the visual band . To express m„ — mo as a function of the angular size, note that flux is proportional 
to 6^T'^, where 9 is the angular size and T is the temperature of the star. 
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Figure 6.1. Spectral density as a function of temperature for several different visual magnitudes 
and observed wavelengths. Atmospheric absorption as a function of the wavelength is not taken 
into account. 
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Figure 6.2. Array configuration used for our analysis. 
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Figure 6.3. A total of 4656 nonzero baselines are available in the array design used in this 
study. Gray scale measures the log of the baseline redundancy. Since the array shown in Figure 
6.2 is almost symmetric with respect to x and y, only a quadrant of the (n, v) plane is displayed. 




Figure 6.4. The four parallel curves indicate the maximum detectable visual magnitude as a 
function of the temperature for several exposure times. Each of these four curves corresponds to 5 
standard deviation measurements and I7I = 0.5. Also shown is the (shaded) region corresponding 
to angular radii between 0.03 mas and 0.6 mas, and observable within 100 hrs of observation 
time. The positions in (T, m^,) space of 2000 stars from the JMMC stellar diameters catalog are 
included. 
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array, we use the JMMC stellar diameters catalog (Lafrasse et al., 2010). We find that ~ 1000 
(out of ~ 33000) stars from the JMMC catalog can be detected within 1 hr, correlations from 
~ 2500 stars can be detected within lOhrs, and ~ 4300 can be detected within 50hrs. In 
Figure 6.2, we show a random sample of 2000 stars (out of 33000) from the JMMC catalog. 
Interstellar reddening may play a role in reducing the number of measurable targets 

6.3 Simulation of realistic data 

Pristine images of disk-like stars, oblate stars, binaries, or featured stars, are first generated. 
The original "pristine" image consists of 2048 x 2048 pixels corresponding to ~ 10 mas x 10 mas 
of angular extension and a wavelength of A = 400 nm. The Fourier transform of the image is 
then performed via an FFT algorithm and normalized so that its value is one at zero baseline. 
This results in a Fourier transform sampled every ~ (8m)/A, i.e., 2 x 10^ cycles per radian 
field-of-view at a wavelength A of 400 nm. We then find the squared-modulus of the degree 
of coherence between the members of each pair of telescopes. This is obtained from a linear 
interpolation of the squared Fourier magnitude in the finely sampled FFT. Diurnal motion is 
not taken into account in the simulations. Diurnal motion plays a significant role in increasing 
the {u, v) coverage when exposure times are long. As a consequence there is less (n, v) coverage 
in the simulations since projected baselines do not drift with time. The smaller {u, v) coverage 
is however compensated by smaller statistical error in the correlation measurements. 

The final step in the simulation phase is the addition of noise to the correlation at each 
baseline. This noise was found to be Gaussian by performing the time integrated product of 
two random streams of simulated photons as detected by a pair of photo-multiplier tubes. The 
standard deviation of the noise added to each pair of telescopes is calculated from eq. 6.2. In 
this study we take the signal bandwidth to be A/ = 200 MHz. An example of simulated data as 
a function of telescope separation is shown in Figure 6.5. This corresponds to a 3'*^ magnitude 
uniform disk star (T = 6000°K) of radius 0.1 mas and lOhrs of observation time. The software 
used for the simulations, as well as the analysis^, was developed in C. 

In section 6.5, the capabilities for reconstructing simple stellar images, with mostly uni- 
form radiance distributions, are discussed in detail. Then the degree of image complexity is 
increased by generating pristine images of stars with nonuniform radiance distributions, e.g., 
limb-darkening and star spots. These simulated images correspond to black-bodies of a specified 
temperature containing an arbitrary number of "star spots" of variable size, temperature, and 
location at the surface of the spherical star in this three-dimensional model. The simulated stellar 



■^See sections 6.4.1 and 6.4.2 for details on the analysis. Some variants of the analysis software were developed 
in MATLAB. All software is available upon request. 
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surface is then projected onto a plane, so that spots located near the edge of the visible half-sphere 

appear more elongated than those located near the center. Additionally, limb-darkening is 
included by assuming that the stellar atmosphere has a constant opacity (more details of the 
simulated images are presented in section 7). Then the image reconstruction capabilities are 
quantified in section 7.1. 



The estimation of derivatives of the Fourier log-magnitude is at the heart of the Cauchy- 
Riemann phase recovery algorithm (section 4.5.2), and is thus greatly simplified when data is 
known on a square grid rather than in a 'randomly' sampled way as is directly available from 
observations. Once simulated data are available (or observations in the future), an analytic 
function is fitted to the data. 

An analytic function can be expressed as a linear combination of basis functions. When data 
f{xi) = |7(xi)p are available at baselines Xi, with uncertainty df{xi), the coefficients of the basis 
functions can be found by minimizing the following 



Each Ofc is the coefficient of a basis function gk- The constant a is a scaling factor, and 
i? is a rotation operator. The scaling factor and rotation angle are found by first performing a 
two-dimensional Gaussian fit. Finding the appropriate scale and rotation angle has the advantage 
of reducing the number of basis elements needed to fit the data. 

Basis functions that tend to zero at very large baselines, where data are scarce (see Figure 
6.3), are ideal. For this reason, we use Hermite functions. There are situations where data are 
more easily fit with a different set of basis functions, e.g., a binary with unresolved members, 
where the data is purely periodic. In such a situation, data do not rapidly tend to zero at large 
baselines, so the Hermite function fit may contain a large number of elements and result in high 
frequency noise where data are scarce.^ The best choice of basis functions may therefore depend 
on the structure of the object. However, for consistency, we use the Hermite fit for all the objects 
that we analyze, and find that it gives reasonably good results. 

The minimization problem can be turned into a linear system by setting the set of partial 
derivatives to zero. We typically start with a small number of basis elements, say eight. 



6.4 Data analysis 
6.4.1 Fitting an analytic function to the data 




(6.3) 



basis set consisting of products of sines and cosines is more suitable in this situation 
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and only increase the number of basis elements if the optimized reduced is greater than some 
predefined value. 

6.4.2 Cauchy-Riemann phase reconstruction 

Since the image is real, the Fourier magnitude is symmetric with respect to the origin in the 
(n, v) plane. Lack of phase information results in the reconstructed image being arbitrary up to 
a translation and reflection. 

6.5 Imaging capabilities for simple stellar objects 

In order to perform a model-independent image reconstruction, the phase of the Fourier 
transform needs to be recovered from magnitude information only (Labeyrie et al., 2006). The 
Cauchy-Riemann phase reconstruction technique is discussed in section 4.5.2, and we will use 
this to obtain the results presented below. 

We investigated the imaging capabilities for simple objects,^ namely uniform disk-like stars, 
oblate rotating stars, binaries, and more complex images. For most of the objects that we 
consider, image reconstruction is not necessary, i.e., from the Fourier magnitude alone, one can 
extract radii, oblateness, relative brightness in binaries, etc. Estimation of these parameters is 
probably more accurate when extracted directly from Fourier magnitude data only, especially if 
some apriori knowledge of the original image is available. However, measuring simple parameters 
from reconstructed images is the first step in quantifying reconstruction capabilities with lACT 
arrays. We assume no apriori knowledge of the images that are being reconstructed, and then 
do a statistical study of the uncertainties of the reconstructed parameters. 

6.5.1 Uniform disks 

In order to quantify the uncertainty in the reconstructed radius, we simulate data corre- 
sponding to 6*'* magnitude stars {T = 6000° K) with disk radii up to 1 mas for 50 hours of 
exposure time^. An example of such a reconstruction is shown in Figure 6.6, where the radiance 
is shown in arbitrary units between and 1. For a uniform disk, the reconstructed phase is null 

in the first lobe, and wc find that the rnis deviations from the true phase arc approximately 
O.lQrad in the null zone. A first look at the reconstruction in Figure 6.6 reveals that the edge 
of the reconstructed disk is not sharp, so a threshold in the radiance was applied for measuring 
the radius. The radius is measured by counting pixels above a threshold and noting that the area 

''For preliminary study see Nunez et al. (2010). 

^ This brightness and exposure time correspond to uncertainties in the simulated data of a few percent. Such 
long exposure times can be reduced to a few hours as is shown in Figure 6.9 
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Figure 6.5. Example of simulated data for a 3^'^ magnitude uniform disk star (T = 6000° /C) of 
radius 0.1 mas and 10 hrs of observation time. Here we show I7I (instead of the directly measured 
I7P) as a function of telescope separation. 
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Figure 6.6. Example of a reconstructed uniform disk of radius 0.1 mas. Also shown is a slice 
of the reconstructed image (solid line) compared to a slice of the pristine image convolved with 
the PSF of the array (dashed line). 
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of the disk is proportional to the number of pixels passing the threshold. After experimenting 
with different radii, we chose the threshold for measuring the radius to be 0.2. We can now 
compare the simulated and reconstructed radii as is shown in Figure 6.7, where each point in the 
Figure corresponds to an individual simulation (including noise) and reconstruction. Further 
optimization in the threshold for measuring the radius should further improve the precision. 

Figure 6.7 clearly shows that stellar radii ranging from 0.03 mas to 0.6 mas can be measured 
with uncertainties ranging between subpercent and a few percent (Figures 6.8 and 6.9). It can 
be seen from Figure 6.7, that the uncertainty increases linearly as a function of the pristine 
(simulated) radius. This is due to a decrease in the number of baselines that measure a high 
degree of correlation when the angular diameter increases. As the pristine radius Q decreases, 
the distance to the first zero in the correlation increases as ^ so the number of telescopes 
contained within the airy disk increases as B~'^. Consequently, decreasing the pristine radius 
is equivalent to increasing the number of independent measurements by a factor of Q"'^. Since 
the uncertainty decreases as the square root of the number of independent measurements, the 
error decreases linearly with the radius. For radii above 0.6 mas, there are simply not enough 
baselines to constrain the Fourier plane information for image reconstruction. For radii greater 
than 0.6 mas, the distance to the first zero in the degree of correlation is of the order of 100 m, 
but only baselines at 35 m and 50 m are capable of measuring the Fourier magnitude with more 
than 3 standard deviations (see eq. 6.2). In Figure 6.9 we show the relative percent error (RMS 
of a statistic) as a function of time for two radii, where it can be seen that a relative error of a 
few percent is achieved after only a few hours. 

6.5.2 Oblate stars 

For oblate stars we use the same magnitude and exposure parameters that are used for 
disk-like stars. Uniform oblate stars can be described by three parameters: the semimajor axis 
a, the semiminor axis 6, and the orientation angle Q. Judging from the limitations obtained from 
reconstructing disks, we produce pristine images whose values for a and h are random numbers 
less than Imas, and choose 1 < a/6 < 2 . The value of the orientation angle Q also varies 
randomly between 0° and 90°. A typical image reconstruction can be seen in Figure 6.10. 

After applying a threshold on pixel values as was done for disk shaped stars, the reconstructed 
parameters are found by calculating the inertia tensor of the reconstructed image. The eigen- 
values and eigenvectors of the inertia tensor provide information for the reconstructed values 
of a, 6 and Q. To do this, we make use of the relation between the matrix eigenvalue and the 
semimajor /minor axes I^x = \o^^i where M is the integrated brightness. A similar relation 
for /„„ holds for the semiminor axis h. 
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Figure 6.7. Simulated vs. Reconstructed radii for magnitude 6 stars with 50 hours of 
observation time (see footnote 5). The top sub-figure shows the uncertainty for a 0.2 mas 
measurement. The bottom sub- figure shows the residual (Reconstructed-Real) along with the 
uncertainty in the radius. 
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Figure 6.8. Histogram of real radius minus reconstructed radius for 50 hours of exposure time 
on a 6*'^ magnitude star (T = 6000 K) of 0.1 mas radius. 
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Figure 6.9. Percent error as a function of time for two reconstructed radii (m^ = 6). This error 
was found by performing several reconstructions for each radius and exposure time, and then 
taking the standard deviation of the reconstructed radius. 
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Figure 6.10. Simulated and reconstructed oblate rotator of magnitude 6 and 50 hours of 
observation time. 
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The resulting reconstructed semimaj or /minor function of their real values have 

a similar structure as the scatter plot for reconstructed radii shown in Figure 6.7, and are 
well reconstructed up to 0.5 mas within a few percent. In Figure 6.11, it can be seen how the 
uncertainty in the reconstructed oblatcness a/h increases with increasing oblateness. As with 
disk shaped stars (section 6.5.1), the uncertainty in the reconstructed semimaj or /minor axes 
decreases as the square root of the number of baselines measuring a high degree of correlation. 
Therefore, the uncertainty in the reconstructed semimajor/minor axes is proportional to y/ab, 
and the error in the reconstructed oblateness is proportional to ~ y^i/F+crY'^- 

The reconstructed orientation angle as a function of the pristine angle is shown in Figure 6.12, 
and several factors play a role in the uncertainty of the reconstructed value. For a fixed value 
of a and b, the orientation of the telescope array with respect to the main lobe of the Fourier 
magnitude determines the number of baselines that measure a high degree of correlation. The 
number of baselines that measure a high degree of correlation is greater when the main lobe of 
the Fourier magnitude is aligned with the x or y direction of the array (see Figure 6.2), and is 
smaller by a factor of ~ \/2 (assuming a uniform grid of telescopes) when its main axis is at 
45° with respect to the array. However, the uncertainty (proportional to spread of points) in 
Figure 6.12 does not appear to be symmetric at 0° and 90°, and is smaller at 90°. This due 
to the way the phase is reconstructed, i.e., due to the slicing of the Fourier plane along the u 
or V directions (see section 4.5.5). In the case of Figure 6.12, the {u,v) plane is sliced along 
the u direction, with a single orthogonal reference slice along the v direction. The main lobe of 
the Fourier magnitude of an oblate star has more slices passing through it when it is elongated 
along the v direction (corresponding to an orientation angle of 90° in image space), yielding a 
better reconstruction. This is in contrast to the orthogonal case of 0°, where the main lobe of 
the Fourier magnitude has a smaller number of slices passing through it. 

6.5.3 Binary stars 

Simulated data arc generated for 5*^^ magnitude binary stars (T = 6000° K), and an exposure 
of 50 hours after noting that the uncertainty in the degree of correlation is of the order of a 
few percent (cq. 6.2). Binaries stars arc parametrized by the radii ri and r2 of each star, their 
separation d, position angle 6, and relative brightness in arbitrary units between and 1. We 
generate pristine images with random parameters within the following ranges: radii are less than 
0.3 mas, angular separations are less than 1.5 mas, the relative brightness per unit area is less 
than or equal to 1, and the orientation angle is less than 90°. A typical reconstruction can be 
seen in Figure 6.13. 

To measure the reconstructed parameters we identify the two brightest spots (regions) whose 
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Figure 6.11. Real vs. reconstructed a/6 for oblate stars. The distance between the two 
linear envelopes is 2 standard deviations. All pristine images that have either a > 0.5 mas or 
6 > 0.4 mas are not included 
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Figure 6.12. Real vs. reconstructed orientation angle for oblate stars. All pristine images 
that have either a > 0.5 mas, b > 0.4 mas or a/b < 1.1 are not included. We also note that 
reconstructed angles are always less than 45° due to the fact that the Fourier magnitude data is 
the same for pristine images flipped about the x, y or x and y axes. Therefore, for all pristine 
angles larger than 45°, we replace the reconstructed angle by ^^g^ = 90° — 9rec- 
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pixel values exceed a threshold of 0.2. We then find the radius for each bright spot and its 
centroid position. Identifying spots is a non trivial task in noisy reconstructed images and our 
analysis sometimes fails to identify the "correct" reconstructed spots. For example, a common 
issue is that close reconstructed spots that arc faintly connected by artifacts, arc sometimes 
interpreted as a single spot. It should be again stressed that image reconstruction may not be 
the best way to measure reconstructed parameters. For example, the data can just as well be 
fit by the general form of the Fourier magnitude of a resolved binary system (containing only a 
few parameters). 

In Figure 6.14, we show reconstructed angular separations as a function of their real values. 
The reconstructed values of the angular separation are found to within ~ 5% of their real values 
and cannot be much larger than what is allowed by the smallest baselines. We find that stars 
separated by more than dmax ~ 0.75 mas are not well reconstructed since the variations in the 
Fourier magnitude start to become comparable to the shortest baseline. 

In Figure 6.15 we show the reconstructed values of the radii as a function of their pristine 
values. We find ~ 10% uncertainties in each of the reconstructed radii. Aside from the angular 
separation, a variable that plays a role in successfully reconstructing pristine radii is the ratio 
of absolute brightness^ of both binary members. When one of the two members is more than 
~ 3 times brighter than the other, the fainter star is found to be smaller than its pristine value, 
and sometimes not found at all when one of the members is more than ~ 10 times brighter 
than the other. This is in part because the sinusoidal variations in the Fourier magnitude 
start to become comparable to the uncertainty. For example: a binary star whose individual 
components cannot be resolved, with one component 20 times brighter that the other, has relative 
variations of ~ 10%. With all the redundant baselines, a few percent uncertainty in the measured 
degree of correlation is sufficient to accurately measure these variations. However, when the 
binary components are resolved, the relative variations decrease with increasing baseline and 
baseline redundancy is not sufficient to reduce the uncertainty in the measurement of the Fourier 
magnitude. This signal to noise issue can of course be improved by increasing exposure time. 

There are also issues related to algorithm performance. One such problem has to do with the 
fit of the data to an analytic function (see section 6.4.1). When the scale of the fit (found by an 
initial Gauss fit) is found to be too small, too many basis elements are used to reconstruct the 
data, and high frequency artifacts appear in reconstructions. Small initial scales are typically 
related to the binary separation as opposed to the size of individual components, and it is 
the latter which correctly sets the scale of the fit. Artifacts may be then mistaken for binary 



^Product of relative brightness per unit area (in arbitrary linear units between and 1) and relative area of 
both stars 
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Figure 6.13. Simulated and reconstructed binary of magnitude 6 and 50 hours of observation 
time. 
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Figure 6.14. Real vs. reconstructed angular separation in binary stars. Binary stars whose 
relative brightness is less than 0.3 are not included in this plot. 
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Figure 6.15. Real versus reconstructed radii for binary stars. We only include cases where one 
of the members is less than 3 times brighter than the other. 
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components, and incorrect reconstructed parameters may be found. Results improve significantly 
when either the correct scale is set or (model-independent) image post-processing is performed 
(see section 7.2). 



CHAPTER 7 



IMAGING STELLAR FEATURES AND 
NONUNIFORM RADIANCE 
DISTRIBUTIONS 

In this chapter, pristine images are generated with varying complexity. In section 7.1, we 
provide two example reconstructions using the Cauchy-Riemann algorithm, and start to quantify 
the reconstruction capabilities. Then we introduce some postprocessing techniques and quantify 
the improvements in a more systematic way: We investigate the reconstruction capabiUties for 
images with increasing complexity by first generating pristine images of stars with limb-darkened 
atmospheres, then we introduce a localized bright or dark feature, and finally increase the number 
of features and explore some of the parameter space, i.e., spot size, location, etc. 

7.1 Featured images with Cauchy-Riemann 
reconstructed images 

We now present two examples of more complex images, and show that the capabilities can, 
to a large extent, be understood from results of less complex images, such as uniform disks and 
binaries. In Figure 7.1 we show the reconstruction of the image of a star with a dark band (an 
obscuring disk for example), corresponding to a 4*^ magnitude star and 10 hrs of observation 
time. The metric used to quantify the agreement with the pristine image (bottom left corner of 
Figure 7.1) is a normalized correlation^ whose absolute value ranges between (no correlation) 
and 1 (perfect correlation/anti-correlation). To quantify the uncertainty in the correlation, we 
perform the noisy simulation and reconstruction several times, and find the standard deviation 
of the degree of correlation. In the case of Figure 7.1, the correlation c is c = 0.947 it 0.001. Note 
that the uncertainty in the correlation is small, and the image reconstruction is not perfect, 
which implies that the reconstruction is not only affected by the SNR level, but also by the 
reconstruction algorithm performance limitations. 



^For two images Aij and Bij, the normalized correlation dj is dj = 
Maxk,i ^N^^{aAO-B)^^^'^j^ (Aij — A){Bi^k.j+i ~ B)^, where cta and ctb are the standard deviations of 
images A and B, and A, B are the image averages. 
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In order to determine the confidence with which we can detect the feature (dark region within 
the disk), we quantify the difference between the reconstruction and a featureless image, i.e., 
a uniform disk whose radius matches the radius of the pristine image. By finding the degree 
of correlation of the reconstruction with a the pristine image (convolved with the PSF), and 
comparing it to the correlation of the reconstruction with a uniform disk, we can quantify the 
confidence level to which our reconstruction is not a uniform disk. We find that the correlation 
of the reconstructed image with a uniform disk is c = 0.880 zb 0.001. This is lower than c = 
0.947 ±0.001 by 61 standard deviations, and this difference is a measure of the confidence with 
which we can establish that the reconstruction does not correspond to that of a featureless star. 

Another example of a complex image reconstruction can be seen in Figure 7.2. Here we 
show the reconstruction of a star with a dark spot, whose correlation with the pristine image 
is c = 0.940 ± 0.001. We can compare this correlation with the correlation of the reconstructed 
image and a uniform disk, which is c = 0.904 it 0.001 (lower by 30 standard deviations). 

For both examples, we also calculate the correlation with the pristine image as a function of 
the angular scale (in mas) of the pristine image. We then find the degree of correlation of each 
reconstruction and its pristine image, and also the degree of correlation of the reconstruction 
with a uniform disk. By comparing these two correlation values, we can estimate the smallest 
feature (spot) that can be reconstructed^. Below some point we no longer distinguish between 
the reconstruction of the featured image and that of a uniform disk. We find that the smallest 
feature that can be reconstructed is close to 0.05 mas. This can already be seen from the order of 
magnitude estimate made in section 6.2 and a comparable value of 0.03 mas is found in section 
6.5.1. When pristine images have angular sizes greater than ~ 0.8 mas, the degree of correlation 
drops significantly due to lack of short baselines. 

7.2 Improved analysis and postprocessing routines 

The resulting reconstructed image with this estimated phase obtained from the Cauchy- 
Riemann algorithm is sometimes not ideal, and so is taken as a first guess for iterative algorithms. 
We have explored the use of the Gerchberg-Saxton (error-reduction) algorithm described in 
section 4.5.7. Recall that constraints must be applied in both the Fourier domain and the image 
domain for this algorithm to converge. The Fourier constraint consists in replacing the Fourier 
magnitude of the image by that given by the data. The constraints in image space can be very 
general. The image constraint that we impose comprises applying a mask, so that only pixels 
within a certain region are allowed to have positive nonzero values. For the images presented 



more careful analysis would require only changing the spot size as opposed to scaling the whole pristine 

image. 
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Figure 7.1. Star with obscuring disk (raw reconstruction). This corresponds to 4*^ magnitude 
and 10 hrs of observation time. The correlation between the real and reconstructed image is 
c = 0.947 lb 0.001. Note that an inverted gray scale is used. 




Figure 7.2. Star with dark spot (raw reconstruction). This corresponds to a 4 magnitude 
star, 10 hrs of observation time and a degree of correlation of 0.940 it 0.001. 
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below, the mask is a circular region whose radius is typically found by measuring the radius of 
the first guess obtained from the Cauchy-Riemann approach. In all reconstructions where the 
Gerchberg-Saxton is used, we perform 50 iterations, and found that more iterations typically do 
not produce significant changes in the reconstruction (Nunez et al., 2012). 

Another postprocessing application that has been utilized is MiRA (Multi-aperture image 
reconstruction algorithm) (Thiebaut, 2009). MiRA has become a standard tool for image 
reconstruction in amplitude (Michelson) interferometry. MiRA is an iterative procedure which 
slightly modifies image pixel values so as to maximize the agreement with the data. In the 
image reconstruction process, additional constraints such as smoothness or compactness can 
be applied simultaneously, but this is something that we have not yet experimented with, i.e., 
the regularization parameter is set to zero for all reconstructions presented here. In the results 
presented below, the number of iterations is set by the default stopping criterion of the optimizer. 
The MiRA software only uses existing data in the (u, v) plane, as opposed to using the fit of an 
analytic function to the data as is done in the Cauchy-Riemann and Gerchberg-Saxton routines. 
This results in removing artifacts in the reconstruction that can be caused by the fit of an 
analytic function to the data. 

A systematic study of the improvements with image postprocessing is presented in sections 
7.2.1, 7.2.2 and 7.2.3. We investigate the performance of each algorithm individually as well 
as the performance of linking algorithms together, particularly the Cauchy-Riemann algorithm, 
followed by the Gerchberg-Saxton and MiRA (Figure 7.3). The order of postprocessing routines 
is also investigated and results are presented in section 7.2.2. 

7.2.1 Limb-dcirkening 

Image reconstruction is actually not necessary for the study of limb-darkening, which can 
be studied directly from the knowledge of the squared degree of coherence. A direct analysis of 
the data, with no image reconstruction, is likely to yield better results than the ones presented 
below. However, it is instructive to first see this effect in reconstructed images before adding 
stellar features to the simulated pristine images. Limb darkening can be approximately modeled 
with a single parameter a as /(0)//o = (cos0)" (Hestroffer, 1997), where is the angle between 
the line of sight and the perpendicular to the stellar surface isopotential. The values of a depend 
on the wavelength and can be found by assuming hydrostatic equilibrium. At a wavelength of 
400 nm, a ^ 0.7 (Hestroffer &; Magnan, 1998) for sun-type stars, and deviations from this value 
may be indicative of stellar mass loss. An example of the reconstruction of a limb darkened star 
with a = 5 is shown in Figure 7.4; such a large value is chosen so that the effect is clearly visible 
in a two dimensional image with linear scale. More realistic values are considered below. To 
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Figure 7.3. Analysis overview for data simulation and image reconstruction. 



obtain Figure 7.4, a first estimate of the phase was obtained from the Cauchy-Riemann algorithm 
and used to generate a raw image. Then the Gerchberg-Saxton postprocessing loop (Figure 4.7) 
was performed several (50) times. In Figure 7.5 the ratio of the average radius at half maximum 
i?i/2 and the nominal radius of Rq, is shown as a function of the limb darkening parameter 
a. Here, data were simulated corresponding to stars with apparent visual magnitude = 3, 
temperature T = 6000° K, radii of Rq = 0.3 mas, lOhrs of observation time and A = 400 nm. 
The ratio Rij2/Ro is less than 1, even in the absense of noise since the reconstruction is at best 
a convolution^ with the array point-spread- function (PSF). 

From Figure 7.5 we can see that we have some sensitivity to changes in the limb-darkening 
parameter a. Stars experiencing high mass loss rates are likely to have high values of a, and if 
we fit a uniform disk function to a limb-darkened reconstruction, the fit yields a smaller radius. 
For example, in the case of a = 2.0, a uniform disk fit yields an angular radius that is smaller 



''a "perfect" reconstruction gives -R1/2/-R0 = 0.96 for a = 0. The extra discrepancy is due to a small hot-spot 
in the reconstruction. R1/2 < Ro since the reconstruction is normalized to the highest pixel value. 
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Figure 7.4. Image reconstruction of a star with limb darkening parameter a = 5, apparent 
visual magnitude = 3 and lOhrs of observation time. The pristine starting image from 
which intensity interferometric data were simulated is shown in the bottom left corner with the 
same contour lines. The Cauchy-Riemann phase reconstruction was performed to produce a raw 
image, and then the Gerchberg-Saxton routine was implemented to produce the postprocessed 
image shown. 

by 7% (still larger than the subpercent uncertainties found in radius measurements (Nufiez 
et al., 2012b)). A real example is the case of the star Deneb, where the difference between the 
extracted uniform disk diameter {9ud = 2.40 it 0.06 mas) and the limb-darkened diameter is 
0.1 mas (Aufdenberg et al., 2002b), and its measured mass loss rate is lO~'^M0yr~^. 

7.2.2 Stars with single features 

Stars were simulated as black bodies with a localized feature of a higher or lower temperature 
as described in section 6.3. In the simulated images, the effect of limb darkening is included as 
described in the previous section. Here the full reconstruction analysis was used, which consists 
in first recovering a raw image from the Cauchy-Riemann algorithm. The raw image is then used 
as a starting point for several iterations of the Gerchberg-Saxton loop (see Figure 4.7), and finally 
the output of the Gerchberg-Saxton algorithm is the starting image for the MiRA algorithm. 
Examples can be seen in Figures 7.6 and 7.7, corresponding to the postprocessed reconstructions 
of bright stars of = 3, lOhrs of observation time and a temperature T = 6000° K. In Figure 
7.6 the temperature of the spot is Tgpot = 6500° K, and in Figure 7.7 the temperature of the 
spot is Tspot = 5500° K. 

We can estimate the smallest temperature contrast that can be detected by varying the 
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Figure 7.5. For each value of a, SII data were simulated corresponding to stars with apparent 
visual magnitude = 3, lOhrs of observation time and A = 400 nm. For each image 
reconstruction, a 1-dimensional profile is found by calculating the average intensity at each 
radial position. The half height i?i/2 (angular radius at half maximum) is found. For each value 
of a, the process was repeated several (10) times, and error bars are the standard deviation of 
the distribution of i2i/2(a)- The ratio R112/R0 < 1 since the image reconstruction is at best a 
convolution with the PSF of the array. 
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Figure 7.6. Reconstructed bright spot. This simulated reconstruction corresponds to a star of 
TTiti = 3, lOhrs of observation time, T = 6000 K, and spot temperature of Tspot = 6500° K. 
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Figure 7.7. Reconstructed dark spot. This simulated reconstruction corresponds to a star of 
= 3, lOhrs of observation time, T = 6000 K, and spot temperature of Tgpot = 5500° K. 

parameters in the model producing the pristine image. The performance in terms of temperature 
contrast obviously depends on several variables such as the size, location and shape of the spot. 
To quantify the smallest detectable spot temperature contrast, we calculate the normalized 
correlation (see footnote 1) between the reconstructed image and the pristine image convolved 
with the array PSF. This correlation is difficult to interpret by itself, so that it is compared with 
the correlation between the reconstruction and a simulated star with no spots. By comparing 
these two values we can have an idea of the confidence level for reconstructing spots with different 
temperatures. This comparison is shown in Figure 7.8, where the top curve corresponds to the 
correlation as a function of spot temperature 6000° K+AT between the reconstructed images and 
the pristine images, and the lower curve is the correlation between the reconstructed images and 
a spotless disk of the same size as the pristine image. A total of 26 stars were simulated, and the 
uncertainty in the correlation was estimated by performing several (10) reconstructions for one 
particular case (AT = 500° K). From the figure it can be seen that spots are accurately imaged 
when AT < —700° K or AT > 200° K approximately. For a black body of spectral density B{T), 
a temperature difference AT < —700° K corresponds to a flux ratio B{T + AT)/ B{T) < 0.45, 
and a temperature difference AT > 200° K corresponds to flux ratios B{T + AT)/B{T) > 1.2. 
This asymmetry can be partly understood in terms of the brightness ratio between black bodies 
B(T + AT)/ B{T), whose rate of change is higher when AT > than when AT < 0. This 
however does not fully account for the asymmetry between cool and hot spots. Most of the 
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Figure 7.8. The top curve and data points correspond to the correlation between reconstructed 
images and their corresponding pristine images containing spots of temperature T + AT, where 
T = 6000° K (my = 3, lOhrs of observation time). The lower curve and data points are 
the correlation between the reconstructed image and a uniform disk of the same size as the 
pristine image. To estimate the uncertainties, we performed several reconstructions and found 
the statistical standard deviation of the correlation. 



asymmetry is due to the fact that all the simulated stars in Figure 7.8 have the same integrated 
brightness, and the radiance per solid angle is larger for a star containing a bright spot than 
for an annular region in a star containing a dark spot. The same analysis can be performed 
by simulating stars with different integrated brightness, but the estimate becomes unnecessarily 
cumbersome and implies knowledge that we would not have access to prior to performing an 
image reconstruction.'^ We should also not forget that this is an estimate, and in a more precise 
calculation we would need to consider additional variables such as spot size, position, etc. 

To test whether the full chain of algorithms is needed to produce Figures 7.6 and 7.7, and to 
investigate algorithm performance, we reconstruct Figures 7.6 and 7.7 with different algorithms 
and combinations of algorithms. Then we calculate the correlation of the reconstructions with 
the pristine image (convolved with the PSF). The results are shown in Table 7.1, and the 
reconstruction for each algorithm combination is shown in Figures 7.9 and 7.10. The correlation 



''For example, we would need to have information on the radiance per solid angle in an annular region in a star 
containing a dark spot 
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Table 7.1. Table comparing correlations of different algorithms. The algorithms are CR 

(Cauchy-Riemann), GS (Gerchberg-Saxton) and MiRA. The first two colums of numbers corre- 
spond to correlations of Figure 7.6, and the second pair of columns corresponds to correlations 
of Figure 7.7. The bold columns correspond to correlations between the reconstructed image 
and the pristine image. The columns not in bold correspond to the correlation with a uniform 
disk from these it can be seen that bright spots are more easily detected. The uncertainty in 
the correlation is Ac = 0.001. 



Algorithm 




C (F 7.6) 




C (F 7.7) 








Pristine 


UD 


Pristine 


UD 


CR 




0.954 


0.928 


0.942 


0.944 


GS 




0.955 


0.928 


0.943 


0.944 


MiRA 




0.978 


0.972 


0.974 


0.980 


CR ^ GS 




0.973 


0.954 


0.966 


0.968 


CR ^ MiRA 




0.979 


0.963 


0.973 


0.971 


GS ^ MiRA 




0.976 


0.965 


0.973 


0.978 


MiRA ^ GS 




0.970 


0.961 


0.965 


0.972 


CR ^ MiRA 


^ GS 


0.968 


0.952 


0.963 


0.961 


CR ^ GS ^ 


MiRA 


0.980 


0.961 


0.977 


0.973 



is found for reconstructions using combinations of Cauchy-Riemann^ (CR), Gerchbcrg-Saxton 
(GS), and MiRA. The single most effective algorithm for these reconstructions is MiRA, and 
the highest correlation is obtained by using the full analysis: Cauchy-Ricmann, followed by 
Gerchberg-Saxton and MiRA. When MiRA or the Gerchbcrg-Saxton algorithms arc used directly, 
the reconstructed image is usually symmetric. The Cauchy-Riemann stage is usually the one 
responsible for roughly reconstructing asymmetries such as the bright or dark spot displayed in 
Figures 7.6 and 7.7. The role of Gerchberg-Saxton is more to improve the phase reconstruction. 
MiRA plays the important role of removing artifacts, caused for example by the data fitting in the 
Cauchy-Riemann phase, and improving overall definition. Even though non-symmetric images 
can be reconstructed, the final product still is somewhat more symmetric than the pristine image 
for reasons that are still under investigation. When the correlation with the pristine image is 
compared to the correlation with a uniform disk, we can again see that the bright spot (Figure 
7.6) is more easily detected than the dark spot (Figure 7.7). 



®It only makes sense to use the Cauchy-Riemann algorithm first, since this is not an iterative algorithm relying 
on a first guess. 




Figure 7.10. Image reconstructions for different algorithm combinations. The pristine corre- 
sponds to that of Figure 7.7. 
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A better estimate of the smallest temperature difference that can be imaged requires an 
exhaustive exploration of parameter space, but temperature differences of less than AT ~ 200° K 
do not seem to be possible to image when the same brightness, temperature, exposure time, 
angular diameter, spot size and spot position as above are used. Results arc likely to improve 
for hotter stars than those simulated above since signal-to-noisc is higher and also the brightness 
contrast is higher for the same relative temperature differences (AT/T). Another question 
related to imaging single spots is that of finding the smallest spot that can be reconstructed. In 
previous work (Chapter 6, (Nuficz ct al., 2012b)), we show that the smallest possible spot that 
can be reconstructed is given by the PSF of the lACT array used in the simulations, namely 
0.06 mas. 

7.2.3 Multiple features 

As a natural extension to the simulations presented above, data were simulated corresponding 
to stars with two or more recognizable features. In Figures 7.11 and 7.12, reconstructions of stars 
containing several hot spots are shown. The brightness and exposure time are the same as those 
used to simulate single-spot stars (m^ = 3, lOhrs). A detailed investigation of reconstruction of 
two-spot stars was not performed, but the general behavior is similar to that presented in the 
section 7.2.2. The reconstructions improve significantly when the pristine image has a higher 
degree of symmetry, e.g., when both spots lie along a line that splits the star in two. This is 
expected since phase reconstruction is not really necessary for centro-symmetric images. For this 
reason, we tested reconstructions with nonsymmetric pristine images. Even though the shape 
of the spots is usually not well reconstructed, the approximate position and size are reasonably 
accurate. 

The reconstruction and identification of features degrades as the number of features in the 
pristine image is increased. A common characteristic of reconstructing stars with several features, 
is that the larger features in the pristine image are better reconstructed. This is more so in the 
case of stars containing darker regions. Nevertheless, information of positions, sizes and relative 
brightness of star spots can still be extracted. In Figure 7.12 a reconstruction of a star containing 
three hot spots of different sizes and relative brightness is shown. This simulated reconstruction 
corresponds to the same brightness and exposure parameters as those of Figure 7.11. 
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Figure 7.11. Reconstructed star with two hot spots. The pristine image has a temperature of 
6000° K, and each hot spot has a temperature of 6500° K. The simulated data corresponds to 
niy = 3 and T = 10 hrs. 
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Figure 7.12. Reconstructed star with three hot spots. The pristine image has a temperature of 
6000° K, and the spots have temperatures of 6500° K (top right spot and left spot) and 6800° K 
(lower spot), The simulated data corresponds to rriy = 3 and T = 10 hrs. 



CHAPTER 8 



EXPERIMENTAL EFFORTS 

There are currently several ongoing efforts that aim to measure correlations of intensity 
fluctuations. There are two main experiments in operation at the University of Utah: A 
laboratory experiment and the StarBase observatory. A short part at the end of this chapter 
describes the StarBase observatory, but most of it describes our laboratory efforts. 

8.1 Laboratory efforts 

The laboratory experiment consists in measuring the angular size of an artificial star using two 
miniature telescopes. The artificial star is a small pinhole (< 1 mm) illuminated by an incoherent 
light source (see sections 8.3.3). The miniature telescopes consist of two photo-multiplier tubes 
(PMTs) are placed 3 m away from the pinhole as shown in Figure 8.1. A beam-splitter is used to 
allow us to effectively place the two detectors at zero baseline. One of the PMTs is movable so 
that the baseline can be changed. The light collecting area of the PMTs is restricted to a couple 
of millimeters in diameter so that the individual PMTs do not resolve^ the pinhole. Light from 
the pinhole travels through a "pipe" and the light detectors are placed inside a metal box^ so 
that they only receive light from the pinhole. 

8.2 Slow control and front end electronics 

The electronics in the laboratory arc set Tip essentially the same way as in the camera used 
for on-sky observations at the StarBase observatory. The camera electronics consist of two 
parts. The slow control electronics provide power^ to the front end, digitize the anode current to 
monitor the DC light intensity (/), provide high voltage to the photo-detector and can be used 
to program parameters of the front end electronics. The Slow control was developed by Jeremy 
Smith, Derrick Kress and Janvida Rou in the University of Utah. The front end electronics 
convert high frequency intensity fluctuations A/ down to the single photon level to analog light 



^Note that when using 400 nm hght, the transverse coherence length is < 10 mm 
^The metal box should also help shield against external signals. 

^AU power sources are external batteries to isolate electronics. 
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pulses which can be transported by optical fiber with minimal bandwidth loss and attenuation 

over great distances. The optical fiber signals are then converted back to electrical signals which 
may be correlated (at the central control). A schematic of the electronics is shown in Figure 8.2. 
The front end electronics were developed for lACTs by Rose et al. (2000), and arc described in 
detail by White et al. (2008). More details of the electronics for SII are given by (LeBohec et al., 



We have experimented with a couple of approaches to measure the intensity correlations: one 
is with an analog system, similar to that used by Hanbury-Brown, and another is with a signal 
digitizer so the correlation is obtained off-line. 



In the analog system, one of the output signals undergoes a periodic polarity flip (phase 
switch) , and then both signals are fed into a analog mixer (multiplier) . The output of the mixer 
is fed into a phase sensitive (lock-in) amplifier, whose reference frequency is provided by the 
signal that controls the phase switching. If there is a correlation between the two signals, then 
the output of the linear mixer displays a periodic change from correlation to anti-correlation at 
the phase switching frequency amplified by the lock-in amplifier. A schematic of this system 
can be seen in Figure 8.3. The time constant of the lock-in amplifier determines the integration 
time of the correlation, and the measurement of |7p can be normalized by finding the DC signal 
provided by the slow control. 

The functionality of the analog system was demonstrated by using a pair of fast blinking light 
emitting diodes (LED), which provided short (20ns) and faint (ImV) light pulses correlated 
between the two channels. The output of the lock-in amplifier as a function of pulse frequency 
is shown in Figure 8.4. The top curve corresponds to the pulses being correlated, and the 
bottom curve corresponds to the signals being anti-correlated (since signals are AC coupled). 
Measurements can be compared to the theoretical prediction, which can be calculated as follows. 
Assume two periodic pulses, with period T, and pulse duration At, defined as 



2010). 



8.3 Correlators 



8.3.1 Analog system 




So when < t < T - At 
si when T - At <t <T. 



(8.1) 



Now, noting that the average signal is 



(s) = So + ^(si - So), 



(8.2) 
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Light source ^ 

Pinhole 




Figure 8.1. Illustration of intensity interferometry laboratory experiment. The baseline can be 
adjusted by moving PMT2. The correlation can be found via an analog system (section 8.3.1) 
or digital system (section 8.3.2). 




Figure 8.2. Outline of electronics. The camera, shown on the left, consists of the PMT, the 
slow control, and the front end electronics. The front ten electronics provide power to the PMT 
and read the anode current from the front end electronics. The front end electronics convert the 
electronic signal to an optical signal, sent to the control station (right), and converted back to 
an electronic signal, where it can be correlated with the signal from another telescope. 
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Figure 8.3. Schematic of analog system for measuring the correlation between input signals Ii 
and l2- The polarity of one of the signals is periodically flipped, so that the output of the linear 
mixer periodically changes sign, and can be detected with a lock-in amplifier. 
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Figure 8.4. Measured correlation for pulses with a duration of 20 ns. The integration time of 
the lock-in amplifier was set to 10 s. The top curve corresponds to the pulses being in phase, 
and the bottom curve corresponds to the pulses being out of phase. Both of these curves are 
compared to the theoretical prediction. 
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then the correlation is 



c = 



(8.3) 




At{T - At) 
7^ ■ 



(8.4) 



(8.5) 



When a time delay AT, equal to the length of the pulse, is introduced between the two 
signals, it is straightforward to show that the (time-delayed) correlation CAt is 



In order to compare the expected behavior of c and cat to measurements, they need to be 
scaled accordingly, i.e., the units of the expected values arc expressed in square volts, whereas 
the read out of the lock-in amplifier is in volts. There is also an offset (of unknown origin) of 
100 fiY in the correlation that needs to be added to the expected curves in order to compare 
them with measurements. The expected curves are shown in Figure 8.4, and deviations from 
the data can be seen when the LED frequency approaches the reciprocal of the pulse width 
((2At)-i = 25 MHz) as expected. 



Recent advancements in high speed data acquisition now allow us to continuously digitize 
data at high frequencies while recording the traces in disk. The correlation between intensity 
signals is then obtained off-line by data analysis. An advantage of this very flexible approach is 
that noise from narrow-band sources, such as cell-phones, motors, etc., can be removed through 
signal processing techniques. Since measuring correlations for thermal light sources requires 
long integration times (S> Is), the disadvantage is that the vast amounts of data generated 
(~ 4Gbs"'^) can cause computational difficulties. Much of the work related to digitizing the 
data efficiently has been done by David Kieda, and is described by (LeBohec et al., 2010). 
The two electronic signals can be sampled at 250 MHz with 12 bit resolution using a National 
Instruments PXIe-5122 high speed digitizer. This digitizing system has proven to be capable of 
streaming data for several hours. 

The functionality of the digitizing system was demonstrated by flashing LED pulses at 1 MHz. 
The signals fed to the LEDs consisted of a floor of light of 820 mV and periodic pulses of 840 mV 
that were 8 ns wide. By averaging many traces recorded with the oscilloscope, we estimated the 
number of additional photons associated with the pulse to be 1/5 per pulse. This means that 



X 2 ' 

CAT = —(>s TP 



(8.6) 



8.3.2 Digital system 
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Figure 8.5. Degree of correlation as a function of the time lag between two digitized LED 
signals. LEDs were flashed at 1 MHz and consist of a floor of light of 820 mV and 8 ns pulses of 
840 mV height. The data sample is 1 s long, and the uncertainty in the correlation was found by 
evaluating the correlation every 4 ms and then calculating the statistical standard deviation. 

every 25 pulses on average we have two photons in coincidence between the two channels. A 
time lag of 200 ns was introduced between the two LEDs so the correlation as a function of the 
time lag should be maximal when signals are brought back in time in the data analysis. This 
maximal correlation at 200ns can clearly be seen in Figure 8.5, which corresponds to Is of data. 
The erratic behavior at short time delays (< 150ns) may be indicative of PMT cross talk or 
external noise. 
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8.3.3 Thermal light source 

The light source must provide incoherent hght with a high spectral density and a manageable^ 
photon count rate. A mercury arc lamp of wavelength 435 nm produces incoherent light similar 
to what can be detected from a hot star with an interferometric filter (AA ^ 10 nm). With such 
a source the coherence time is clearly much shorter than the electronic resolution time, and in 
this regime, the SNR given by equation 3.39, i.e., 

SNR = NAa\J\^^/ToAf/2. 

To estimate the SNR, we can first find the photon rate NAaSu by reading the typical anode 
current^ from the PMT photo-cathode, which is typically of the order order of / 10 /xA. The 
gain of the PMTs is of the order oi G ^ 10^ , therefore, the number of photons per unit time is 

NAaAu « « 6 X 10^ s-^ (8.7) 

Since the optical bandwidth is Au = cAA/A^ lO^^Hz and the electronic bandwidth is 
(pessimistically) A/ 100 MHz, then a SNR Ri 3 requires an integration time of ~ 3min. 

Several attempts have been made to measure correlations with a thermal source. In these 
attempts, the analog system was used as described in section 8.3.1. Measuring intensity correla- 
tions has proven to be difficult with a thermal source. It is advantageous to perform an intensity 
interferometry experiment where long integration times are not needed. The electronic time 
resolution or the spectral density cannot be significantly improved. However, in the next section 
we describe an experiment in which an artificial light source with an extremely long coherence 
time is used for intensity interferometry experiments. 

8.3.4 Pseudo-thermal source 

An incoherent light source with a coherence time which is much longer than the electronic 
resolution has been constructed as proposed by Martienssen Sz Spiller (1964). A pseudo-thermal 
light source can be easily created by scattering coherent light (e.g., laser light) off a medium 
that continuously changes with time, therefore producing a time-varying speckle pattern. This 

can be accomplished by shinning laser light (A = 534 nm) through a rotating sheet of ground 
glass as shown in Figure 8.6. The rotating ground glass is put as close as possible to the pinhole 

*A very low count rate would require very long integration times to detect a correlation (see SNR estimate 
below). 

^The slow control actually provides a DC voltage which can be converted to a current by finding an equivalent 
resistance of the electronics. Using a signal generator to provide a calibration signal, we have found the equivalent 
resistance to be 12.2 kfl. 
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to create the artificial star, and a rapidly changing speckle pattern can be seen where the light 
is detected by the PMTs. 

The typical time over which the speckle pattern changes sets the coherence time of the 
source, and can be adjusted by changing the speed of the ground glass, or by using a different 
sized glass grain. The photon detection statistics generated by such a source are very similar 
to those of a thermal light source (Morgan &; Mandel, 1966), except that the detected light 
fluctuation are now dominated by the wave noise instead of the shot noise (see section 3.1.2), 
that is, the fluctuations are proportional to the number of detected photons per resolution time. 
The coherence area is essentially the size of a speckle, and when photo-detectors are separated 
by less than the typical size of a speckle, i.e., within the coherence volume, then their signals are 
correlated. A pseudo-thermal source then allows us to do an intensity inter ferometry experiment 
in "slow motion." 

8.4 Results with pseudo- thermal light source 

Correlations were measured between two PMTs receiving light from an artificial star (pinhole) 
that is illuminated by pseudo-thermal light. 

8.4.1 Individual signals from each channel 

A sample trace is shown in Figure 8.7. The top curve in Figure 8.7 is the raw data obtained 
by the digitizer. Since speckles have a typical duration of a fraction of a millisecond (described 
below in section 8.4.2), the measurement of the normalized degree of correlation is maximal 
when the envelopes of the traces are used as opposed to the raw traces. Since the measured 
signals are AC coupled, the absolute value is taken (middle curve in Figure 8.7) before a low 
pass filter is applied to find the envelope. A filtered trace envelope is displayed in the bottom of 
Figure 8.7. The filtering time constant must be large enough so as to maximize the correlation, 
but still much smaller than the coherence time. Therefore, it is useful to measure the coherence 
time before finding the normalized degree of correlation. 

8.4.2 Temporal coherence 

In order to find the coherence time of the pseudo-thermal light source, the correlation is found 
as a function of a time displacement introduced between the two channels (at zero baseline). 
Since data were taken with the digital system described in section 8.3.2, the time delay is 
introduced off-line. Correlations corresponding to a small sample of 2nis arc shown in Figure 
8.8, where the temporal coherence is measured for both the pseudo-thermal light source, and 
the laser, i.e., rotating ground glass compared to nonrotating ground glass. 
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Figure 8.6. Laser scattered light tiirougii rotating ground glass produces a pseudo-thermal 
light source. 
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Figure 8.7. The top two traces correspond to the raw data obtained with the digitizing system 
in two separate channels. The signal is plotted every 1000 points so that the plot is not saturated. 
In order to maximize the degree of correlation, the envelope of each curve is found before the 
correlation is calculated. To find the envelope, the absolute value of the signal is found before a 
lo pass filter (r < 0.01 ms) is applied. The bottom figure shows the two filtered signals, which 
are clearly correlated. 
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An ideal laser should not display temporal correlations since it has negligible intensity 
fluctuations and emits Poisson distributed light (section 3.1). The laser used in this experiment 
actually contains several modes, with a mode spacing of 303 MHz, and should experience intensity 
modulations with this frequency, but these arc under-sampled at 250 MHz with the digitizing 
system used in this experiment. Even though an ideal laser emits coherent light, intensity 
correlations (or their absence) do not have a relation to the angular radiance distribution of the 
source.^ 

In Figure 8.8, the "photon-bunching" region is clearly seen for delay times smaller than 
^ 0.4 ms, and the coherence time r can be estimated by fitting a Lorentzian, which yields 
r = 0.23 ms. Therefore, a filtering time constant of 0.01 ms is appropriate in order to maximize 
the degree of correlation (see section 8.4.3 for more details). On a thermal source, the type of 
plot shown in Figure 8.8 contains information on the spectrum, and the inverse of the coherence 
time would essentially be the optical bandwidth of the source. For the case of the pseudo- 
thermal source, this contains information on the time-scale in which the speckles pass through 
the detector. If the speckles only originate from the pinhole, then they should all be of similar 
size, i.e., of the size of the coherence area, and therefore, the curve shown in Figure 8.8 should 
have a slope of zero near the origin. The fact that the curve does not level out implies that we may 
also be observing smaller speckles. Smaller speckles would need to originate from sources that 
have a large angular size, and therefore it is possible that they originate from internal refiections 
in our experimental set up (e.g., the metal box where the PMTs are placed). However, we have 
not established this rigorously. 

8.4.3 Uncertainty in the correlation 

A 4000 ns sample trace is shown in Figure 8.9, where individual photons can be seen. At these 
short time-scales no correlation is expected, but by counting the number of photons per unit time 
((n) f« 8 X lO^s-i), the expected SNR in a correlation measurement (with the pseudo-thermal 
source) can be estimated as follows. Recall from section 3.1.2 that the variance of the number 
of photo-clcctrons n contains a Poisson term (shot noise) and a term that is proportional to the 
variance of the time integrated light intensity fj, (eq. 3.20), i.e., 

(An^) = (n) + (Afi^) , 

where 



®The relation between intensity correlations and the degree of coherence (eq. 3.33) only exists in light whose 
electric field tends towards being a Gaussian random variate. 
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Figure 8.8. Temporal correlation for pseudo-thermal light compared to laser light. A time delay 
between the two (nonfiltered) PMT signals was introduced off-line. The "photon-bunching" 
region can be seen below ~ 0.6 ms. The estimation of the error-bars is described in section 8.4.3. 




-80 I ' ' ' ' ' ' ' 1 

500 1000 1500 2000 2500 3000 3500 4000 

time (ns) 

Figure 8.9. Digital counts for each channel. According to this trace, the number of pho- 
tons per unit time is ~ 8 x 10'' s~^. Since the electronic resolution time is 4 ns, there are 
8 X 10^ s~^ X 4 ns = 0.3 photons per resolution time. No obvious correlation can be seen a these 
time-scales. 
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H = ri 




I{r,t')dt' 



In this expression, / is the light intensity integrated over the optical bandwidth and 77 is a 
constant that characterizes the detector. When the electronic time resolution T is much smaller 
than the coherence time, the integrated light intensity is approximately n ~ r]I{t)T. Since there 
is essentially no light between adjacent speckles, the variance of the integrated light intensity is 
approximately 



Therefore, the SNR achieved in one resolution time, is equal to the number of photons per 
resolution time, i.e., rj {I{t)) T = NAaAuT ^ 0.3 from Figure 8.9. For an observation time Tq, 
the SNR increases as the square root of the observation time, i.e.. 



Therefore, when we integrate for ~ Is, we should expect a SNR ~ 4 x 10^ . This value 
should be compared to the measured SNR, described below. 

In order to evaluate the uncertainty in the correlation, the data run is subdivided in many 
small time windows. In each time window, the correlation is found, and then the statistical 
standard deviation of the correlation is found. The duration of each time subdivision must 
be much longer than the coherence time so that it includes the passing of many speckles. To 
find the optimal duration of each time subdivision, we perform a study of the ratio of the 
degree of coherence and its standard deviation {SNR = C/AC) as a function of the number 
of sub-divisions (Figures 8.10 and 8.11) . In Figure 8.10, we can see that the SNR is of the 
order of ~ 30 when the time averaging time window is much longer than the coherence time and 
when more than a few time samples are used (for sufficient statistics). Since the coherence time 
is 0.2 ms, a time window of 2 ms is appropriate, and yields a SNR of ~ 30. The fact that the 
measured SNR is two orders of magnitude smaller than expected indicates that this experiment 
is limited by electronic noise rather than photon noise. 

When no filtering is applied to the data (red curve in Figure 8.10), the SNR decreases as 
soon as the time window is smaller than the coherence time. It is surprising to see that when 
no filtering is applied (red curve in Figure 8.10), the signal, as well as the SNR, increase when 
the length of the time window approaches the electronic resolution (Figure 8.11). The electronic 
resolution time is actually close to the laser mode spacing of 303 MHz, which is under-sampled 
at 250 MHz with the digitizing system. These high frequency correlations are possibly due to 




SNR = NAaAuAT\^f^/Th/T. 



(8.9) 
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Figure 8.10. Ratio of the correlation its statistical standard deviation as a function of the 
number of time subdivisions (oc 1/Time window) in a 0.2 s run. The red curve corresponds to 
the the SNR with no filtering applied to the data. The blue curve has been obtained by applying 
a low-pass filter with r = 0.01 ms. The purple curve corresponds to no filtering, but a time delay 
of 8 /US was applied to test for possible cross talk between PMTs. 
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beating between modes. The correlation of pure laser light (nonrotating ground glass) displays 
the same behavior at very small time windows (light blue curve in Figure 8.11). However, we 
have not rigorously proven this hypothesis and we cannot rule out high frequency noise pick-up 
in the electronics. 

8.4.4 Simulation of slow electronics 

We also study the behavior of the SNR (C/ AC) as a function of the filtering time constant. 
Recall that when the coherence time is much shorter than the electronic resolution time, the 
SNR is diluted by the ratio of the coherence time and the electronic resolution time. When 
the electronic time resolution is much smaller than the coherence time (current regime with 
the pseudo-thermal light source), the SNR is independent of the electronic time resolution. By 
increasing the filtering time constant, the electronics can be made artificially slower. 

Figure 8.12 shows the behavior of the SNR as a function of the filtering time constant. When 
the filtering time constant is much smaller than the coherence time, the SNR is essentially 
constant. When the filtering time constant is larger than the coherence time, the SNR is inversely 
proportional to the Filtering time constant as expected. This is experimental evidence of what is 
described in section 3.3, which essentially says that when electronics are "slow" , or light intensity 
fluctuations are very fast (such as with a thermal source) , then the ratio of the wave-noise (signal) 
and the shot noise is proportional to the ratio of the coherence time and the electronic resolution 
time. 

8.4.5 Spatial coherence and angular diameter measurements 

The most relevant result for SII is the measurement of spatial coherence. The diameters 
of different pinholes emitting pseudo-thermal light were found by measuring the normalized 
intensity correlation as a function of the PMT separation. Since the coherence time is much 
longer than the resolution time, we in fact find the correlation between the envelopes of the 
signals. To find the envelope, take the absolute value and then filter out frequencies larger than 
25 kHz from the individual signals. 

It is reasonable to assume that the pinholes are uniformly illuminated. The diameter of 
a uniformly illuminated pinhole can be extracted by fitting the correlations to an Airy function, 
and by recalling that A9 = 1.22X/D, where D is the first zero of the Airy function. The pinhole 
diameters have nominal values of 0.2 mm, 0.3 mm and 0.5 mm, although the 0.5 mm pinhole was 
made by piercing a sheet of aluminum foil with a metal wire and is therefore only approximately 
this size. Data are presented in Figure 8.13 along with the best fit curves, and the measured 
pinhole diameters are 0.17 ± 0.02 mm, 0.24 it 0.03 mm and 0.32 it 0.05 mm. 
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Figure 8.11. Correlation and its statistical standard deviation as a function of the number 
of time sub-divisions (oc 1/Time window) in a 0.2 s run. The red curve corresponds to the the 
SNR with no filtering applied to the data. To obtain the purple curve, a time delay of 80 ns was 
applied to test for possible cross talk between PMTs. Also shown is the correlation for pure laser 
light (nonrotating ground glass), which displays the same behavior at small time windows as the 
pseudo-thermal light, and shows no significant correlations when the averaging time window is 
much longer than the electronic resolution time. 
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Figure 8.12. Ratio of the degree of correlation and its standard deviation as a function of the 
filtering time constant. The averaging time window is 8 ms. 
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8.5 Image reconstruction from autocorrelation data 

In the experiment with the pseudo-thermal hght souree, time varying speekle patterns can 
be seen at the detectors, and correlations can be measured when both detectors lie on average 
within the same speekle. In view of the experiment with the pseudo-thermal light source, we 
can understand a thermal light source as a completely rough diffracting screen at each instant 
in time. 

In the case of the pseudo-thermal light source experiment, the time integrated correlation (as 
a function of PMT separation) is essentially an autocorrelation of the speckle pattern. Therefore, 
one can in principle take an image of a speckle pattern, compute its autocorrelation, which is 
equivalent to the squared modulus of the Fourier transform of the artificial star, and apply the 
methods developed in Chapter 4 for reconstructing the image of the artificial star. 

This type of experiment is currently in progress as shown schematically in Figure 8.14. An 
artificial star is created by shinning laser light through a mask of any desired shape. Then the 
spatial coherence is broken by scattering the light off a rough (paper) screen. The resulting 
speckle pattern does not change with time, and can then be recorded with a CCD camera. 

Some preliminary results can be seen in Figure 8.15. In this example, a triangular-shaped 
mask is used. The figure shows the speckle pattern, its autocorrelation and corresponding image 
reconstruction. To obtain this image, the Cauchy-Riemann algorithm was applied, followed 
by 100 iterations of the Gerchberg-Saxton algorithm'^. The experiment, which has been done 
in collaboration with Ryan Price and Erik Johnson, will be performed more carefully so that 
actual angular scales of the image can be obtained. 

8.6 Starbase 

As a first test toward implementing SII with lACT arrays, pairs of 12 m telescopes in the 
VERITAS array at the Fred Lawrence Whipple Observatory in Arizona were interconnected 
through digital correlators (Dravins &; LeBohec, 2008). These tests were made during nights 
shared with VHE observations with a very temporary setup and established the need for a 
dedicated test bench on which various options of secondary optics and electronics could be 
evaluated in a realistic environment. In order to satisfy this requirement, the two StarBase 
telescopes were deployed on the site of the Bonneville Seabase diving resort in Grantsville, Utah, 
40 miles west from Salt Lake City. The two telescopes (Figure 8.16) are on a 23 m East- West base 
line. The telescopes had earlier been used in the Telescope Array experimentl operated until 
1998 on the Dugway proving range. Each telescope is a 3 m, f/1 Davies-Cotton light collector 



The constraint applied in image space is that in each iteration, all pixels below 0.02/1.0 are set to zero. 
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Figure 8.13. Normalized degree of correlation as a function of PMT separation for three 
different pinhole sizes. The curves corresponding are the best fit curves. 



Mask 



jaser 



^ea^^xgand^l 



Rough surface 
1 1 1 la ■ ■ ■ f 



Figure 8.14. The laser beam passes through a beam expander and then through a mask of any 
desired shape (e.g., single or multiple circular holes). The rough paper screen reflects a speckle 
pattern which is captured with a CCD camera. The autocorrelation of the speckle pattern is 
the squared magnitude of the Fourier transform of the image on the rough screen. 
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Figure 8.15. The autocorrelation is found for a speckle image obtained with a triangular-shaped 
mask. Image reconstruction techniques allow us to obtain the righ-most image. 
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composed of 19 hexagonal mirror facets ~ 60 cm across. This design is typically used for lACT 
and secondary optics tested on the StarBase telescopes could be used directly on the VERITAS 
telescopes for larger scale tests. 

8.6.1 Sensitivity and outlook 

Using Equation 3.39 with conservative parameters for the StarBase telescopes {A = 6m^ , 
a = 0.2 and A/ = 100 MHz), the 5 standard deviation measurement of a degree of coherence 
|7((i)p = 0.5 will require an observation time T ~ 10minx2.5^™" where rriy represents the visual 
magnitude and where we have made the crude approximation n = 5 x 10^ x 2.5~™"m^s^Hz^ . 
This corresponds to one hour for = 1 and 6.5 hours for = 2 and when considering the 
measurement of |7(d)p ~ 1, these observation times should be divided by four. 

The first objective will be the detection of the degree of coherence for an unresolved object 
(|7(o?)p ^ 1). The distance between the two telescopes being 23m (smaller baselines can be 
obtained during observations to the east and to the west due to the projection effect), at A = 
400 nm the stars have to be smaller than typically ~ 3 mas in diameter. An essentially unresolved 
star suitable for calibration should be less than ~ Imas in diameter. Good candidates for 
this are in increasing order of magnitude a-Leo, 7-Ori, /3-Tau or even ?7-UMa which, should 
be observable as a unresolved object for calibration within 50 minutes LeBohec et al. (2010). 
Alternatively, it will be possible to measure any star as an unresolved object by correlating the 
signals from two channels mounted on the same telescope by means of the camera beam splitter. 
These observations should allow us to establish methods for adjusting the signal time delays 
optimally and also to identify the most effective correlator. A next phase will be dedicated 
to the measurement of a few bright stars in order to further demonstrate the technique. This 
second phase will possibly include the observation of coherence modulation resulting from orbital 
motion in the binary star Spica with a 1.5 mas semi major axis and = 1.0, or even, possibly 
Algol (2.18 mas semi major axis, m^, = 2.1). 

8.6.2 First data sample 

The StarBase observatory recently began taking data of the binary star Spica. Figure 8.17 
shows the digitized signal of each telescope. We have so far only analyzed 1 s of data and 
calculated the normalized correlation as a function of a time lag between the two signals as 
shown in Figure 8.18. At this very preliminary stage we do not expect to sec a signal since an 
integration time of a few hours is necessary to detect a correlation (eq. 3.39). From Figure 8.18, 
we note that the fluctuations in the degree of correlation are comparable to what was obtained 
in the laboratory with the LED experiment when the LEDs were out of phase (section 8.3.2). 
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Figure 8.16. The StarBase 3 m telescopes are protected by buildings which can be rolled open 
for observation. The control room is located in a smaller building located between the two 
telescopes. 
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Figure 8.17. The first individual signals obtained from the binary star Spica. 
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Figure 8.18. Degree of correlation as a function of the time lag between the two channels. This 
corresponds to 1 s of data. 
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A next test at StarBase could be to perform the same LED experiment as was done in the 

laboratory. 



CHAPTER 9 
CONCLUSIONS 

Gamma-ray astronomy is typically not used for stellar astrophysics, but can still further our 
understanding of some isolated stellar systems. This is the case of LSI + 61°303, consisting 
of a main sequence Be star and a compact object, which has been detected in the TeV range 
with VERITAS. This object showed a clear intensity modulation as a function of the orbital 
phase. We describe a gamma-ray attenuation model and apply it to this system. With this 
model, we are able to constrain fundamental parameters of the system such as the mass of the 
compact object and the density of circumstellar matter around the Be star. However, important 
details of this source, such as the circumstellar matter distribution, may only be obtained from 
model-independent imaging at high angular resolution. Interferometry observations have already 
allowed us to reconstruct high angular resolution images of a few stars at optical wavelengths, 
but most bright and nearby stars cannot be resolved with current facilities. 

The recent success of gamma-ray astronomy, as well as the advancements in instrumentation 
and computing technology since the days of the Narrabri intensity interferometer, have prompted 
a revival of optical SII. Kilometric scale arrays of many large light collectors will allow an 
improvement in angular resolution by nearly an order of magnitude when compared with current 
optical amplitude interferometers. Several thousand stars will be resolved with such large arrays 
used as intensity interferometers, and due to the very dense sampling of the (u, v) plane, new 
mathematical (phase retrieval) algorithms will allow for high angular resolution images to be 
reconstructed from SII data. 

We performed a simulation study of the imaging capabilities at 400 nm of an lACT array 
consisting of 97 telescopes separated up to 1.4 km. This is a preliminary design for the CTA 
project, expected to be operational in 2018. Our method uses a model-independent algorithm 
to recover the phase from intensity interferometric data. We tested the method on images 
of increasing degrees of complexity, parameterizing the pristine image, and comparing the 
reconstructed parameter with the pristine parameter. We now summarize our results and briefly 
discuss how fundamental stellar parameters can be constrained with the methods described in 
this thesis. 

We found that for bright disk-like stars (m„ = 6, T > 6000° K), radii are well reconstructed 
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from 0.03 mas to 0.6 mas. Even though using a phase retrieval approach to recover images 
might not be the most efficient way to measure stellar radii, such a study starts to quantify 
the abilities of measuring other scale parameters in more complicated images (e.g. oblateness, 
distance between binary components, star spots, etc.). The range of angular radii that can 
be measured with a CTA-like array (0.03 — 0.6 mas) will complement existing measurements 
(2 — 50 mas) (Haniff, 2001). With the aid of photometry, the effective temperature^ scale of 
stars within 0.03 — 0.6 mas can be extended. 

Binary stars are well reconstructed when one of the members is not much brighter (three 
times as bright) than the other, and when they are not too far apart (< 0.75 mas). As with 
amplitude interferometry, SII, along with spectroscopy, will allow us to determine the masses and 
orbital parameters in binary stars. If measured with enough precision (< 2%) (Andersen, 1991), 
the determination of the mass can be used to test main sequence stellar models. An advantage 
of using an array such as the one used in this study, is that individual radii can be resolved. An 
interesting phenomena to be studied with interacting binary stars is mass transfer (e.g., Zhao 
et al. (2008)), and capabilities for imaging this phenomena can be further investigated. 

For oblate stars, results similar to those obtained for disk-like stars are found. Due to 
the relative ease of SII to observe at short (~ 400 nm) wavelengths, measuring fundamental 
parameters of hot B type stars is possible. B stars are particularly interesting since rapid 
rotation, oblateness, and mass loss are a common feature. We show that oblateness can be 
accurately measured, and the next step is to quantify the capabilities of imaging realistic surface 
brightness distributions in hot stars. By imaging brightness distributions, we will be able to 
study effects such as limb darkening and mass loss in hot massive stars (Ridgway et al., 2009), 
as well as gravity darkening (von Zeipcl, 1924). 

To image nonuniform brightness distributions wc recur to postprocessing routines that signif- 
icantly increase the quality of the reconstructed images. Since the Cauchy-Riemann algorithm 
provides a reasonable first guess of a reconstruction, postprocessing routines consist of convergent 
iterative algorithms such as the Gerchberg-Saxton algorithm and the MiRA (Multi-aperture 
image Reconstruction Algorithm). The postprocessing significantly improves the image recon- 
struction, but the postprocessing routines by themselves are usually not sufficient for performing 
reconstructions, especially when the pristine image is not centro-symmetric. 

A study of the imageability of limb-darkening, which is related to the mass loss rate in hot 
stars (r ~ 10,000° K), indicates that realistic mass loss rates of the order of lO~^M0yr~"^ can 
be imaged. By simulating data corresponding to stars containing bright or cool spots, we find 

^Defined as Te// = { iJ^i^ Y^'^ ^ where R is the radius, L is the luminosity, and a is the Stephan-Boltzmann 
constant. 
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that features with sizes of ^ 0.05 mas and temperature differences of a few hundred °K can be 

reconstructed. Photometric variabihty (at the ~ 10% level) that coincides with the rotational 
period would provide indirect evidence of the presence of such stellar hot spots, and in fact, this 
type of photometric variability has been reported in the past (Balona ct al., 1987). The origin of 
such variability is still highly debated, and actual images would provide another mean to study 
this phenomenon. 

Several experiments are currently performing intensity correlation measurements as a prepara- 
tory stage for the use of lACT arrays as intensity interferometers. These efforts include the 
StarBase observatory and a laboratory experiment in the University of Utah. The Starbase 
telescopes have recently begun taking data and testing electronics. An intensity interferome- 
try experiment that measures correlations from simulated stars emitting pseudo-thermal light 
(essentially incoherent light with long coherence times compared with thermal light) has been 
performed. This experiment allowed us to measure angular diameters as well as to see the 
limitations of the electronics. However, the value of this experiment lies in that it enables an 
intuitive understanding of intensity interferometry in terms of speckles, i.e., when observing 
stars at a narrow optical bandwidth, detectors that are closer to each other than the typical 
size of a speckle, display a higher degree of correlation since they are on average detecting light 
corresponding to the same speckle. 

There are opportunities to further investigate SII simulations and experiments. For example, 
the simulated data used in this study to reconstruct images assumes point-like telescopes and 
does not include the effects of electronic noise or night-sky background. Some of these effects 
arc currently being investigated (Rou, 2012). It is also possible to use more realistic pristine 
images that may be provided from Montc-Carlo radiative transfer models, such as those used 
for modeling the atmospheres of hot Be stars (Carciofi & Bjorkman, 2006). Reconstructions 
from these pristine images may be used to investigate how radiative transfer models can be 
constrained. On the experimental side, we would like to make further attempts to measure 
correlations with an actual thermal source, using both the analog and digital systems. The 
StarBase telescopes have recently begun operating, and correlation measurements from small 
stars and binary systems will soon be performed. The lessons learned from these experiences 
will allow us to achieve the ultimate goal: To view stars as not mere unresolved point sources, 
but as the fascinating objects they truly are. 
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Thesis summary 

Towards optical intensity interferometry for high angular resolution stellar 

astrophysics 

cn . . 

1^ , Imaging Air Cherenkov Telescope (lACT) arrays are primarily used for high energy (TeV) gamma-ray 

■ astronomy, and to study some of the most violent phenomena in the universe. These arrays can also 

be used for doing stellar astrophysics, in some cases by directly analyzing gamma-ray data, but also by 
using them as intensity interferometers at optical wavelengths. This study discusses both approaches 
but focuses on Stellar Intensity Interferometry (SII) virith lACT arrays. 



The binary LS I + 61°303, which consists of a hot Be star and a compact object whose nature is 
i-S^ ■ still subject of debate, has been observed with the Very Energetic Radiation Telescope Array System 

(VERITAS) lACT 1 . The data has been used in the context of a simplified model accounting for the 
modulation of the TeV light-curve [2] in terms of gamma-ray attenuation mechanisms in the circum- 
stellar environment of the Be member of the pair. This has allowed to constrain quantities such as 
the mass, density of circumstellar material, and orbital parameters. Results suggest that the compact 
companion is a black hole as opposed to a neutron star. However, this analysis finds a high column 
density of circumstellar material which conflicts with X-ray observations, suggesting that gamma-ray 
attenuation may not play an important role in the intensity modulation. 

Many obstacles in understanding stars can be overcome with high angular resolution techniques. 
Model independent measurements of radii, scale features, and binary orbital parameters are still needed 
to provide effective constraints on stellar evolutionary models. High angular resolution imaging of fast 
rotating stars (e.g. Be stars) will allow to study effects such as gravity darkening and mass loss. Map- 
ping brightness distributions across the stellar surface will allow to quantify radiatively driven mass 
loss, and imaging interacting binary systems will allow to quantify mass transfer. Direct imaging of 
many of these effects is highly relevant for understanding systems such as LS I + 61°303. 

Recent proposals have been advanced to apply Imaging Air Cherenkov Telescope (lACT) arrays to 
optical Stellar Intensity Interferometry (SII) [Sj. SII with lACT arrays, composed of nearly 100 tele- 
scopes, will provide means to measure fundamental stellar parameters and also open the possibility of 
model-independent imaging. A main limitation of image recovery in intensity interferometry is the loss 
of phase of the complex degree of coherence during the measurement process. Recent developments of 
phase retrieval techniques are used in this study to recover images from simulated data corresponding 
to bright stars {m^ ~ 6) and exposure times of a few hours [HIS]. Scale features such as diameters, 
oblateness, and overall shapes, are reconstructed with uncertainties of a few percent . The capabilities 
for reconstructing stellar images with non-uniform radiance distributions are also quantified. Simula- 
tions show that hot stars (T > 6000° K) containing hot and/or cool localized regions (AT ^ 500° K) 
as small as ^ 0.1 mas can be imaged at short wavelengths (A = 400 nm) [B]. A discussion of current 
experimental efforts for measuring intensity correlations in a laboratory set-up as well as with a pair 
of 3 m telescopes separated by 23 m is presented. 
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Imaging Air Cherenkov Telescope (lACT) arrays are primarily used for high energy (TeV) gamma-ray 
astronomy, and to study some of the most violent phenomena in the universe. These arrays can also be 
^ J used for doing stellar astrophysics, in some cases by directly analyzing gamma-ray data, but mostly by 

using them as intensity interferometers. This study discusses both approaches but focuses on Stellar 
Intensity Interferometry (SII) with lACT arrays. 

The first approach of analyzing TeV data has been used for studying high energy emitting binary 
systems in the context of a simplified model, and has resulted in a method to constrain macroscopic 
parameters of these systems [1] . This method has been applied to the particular case of the binary 
LS I -|-61°303 in [1], which consists of a hot Be star and a compact object whose nature is still subject 
of debate. This source has seen a clear TeV intensity modulation as a function of the orbital phase 
[1] , and an analysis of the TeV light-curve in the context of a model of gamma-ray attenuation mech- 
anisms has allowed to constrain quantities such as the mass, density of circumstellar material, and 
orbital parameters. Results suggest that the source is a black hole as opposed to a neutron star pj. 
However, this analysis finds a high column density of circumstellar material which conflicts with X-ray 
observations, suggesting that gamma-ray attenuation may not play an important role in the intensity 
modulation. 

Many obstacles in understanding stars can be overcome with high angular resolution imaging. Re- 
cent proposals have been advanced to apply Imaging Air Cherenkov Telescope (lACT) arrays to optical 
Stellar Intensity Interferometry (SII) [3]. SII with lACT arrays, composed of nearly 100 telescopes, 
will provide means to measure fundamental stellar parameters and also open the possibility of model- 
independent imaging. A main limitation of image recovery in intensity interferometry is the loss of 
phase of the complex degree of coherence during the measurement process. Recent developments of 
phase retrieval techniques are used in this study to recover images from simulated data corresponding 
to bright stars (m„ ~ 6) and exposure times of a few hours [11 [5]. Scale features such as diameters, 
oblateness, and overall shapes, are reconstructed with uncertainties of a few percent . The capabilities 
for reconstructing stellar images with non-uniform radiance distributions are also quantified. Simula- 
tions show that hot stars (T > 6000° K) containing hot and/or cool localized regions (AT ~ 500° K) 
as small as ^ 0.1 mas can be imaged at short wavelengths (A = 400 nm) [6]. A discussion of current 
experimental efforts for measuring intensity correlations in a laboratory set-up as well as with a pair 
of 3 m telescopes separated by 23 m is presented. 
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Optical stellar intensity interferometry with air Cherenkov telescope arrays, 
composed of nearly 100 telescopes, will provide means to measure fundamental 
stellar parameters and also open the possibility of model-independent imaging. 
In addition to sensitivity issues, a main limitation of image recovery in intensity 
interferometry is the loss of phase of the complex degree of coherence during 
the measurement process. Here we make use of recent developments of phase 
retrieval techniques to recover images from simulated data. For bright stars 
(niy ^ 6) and exposure times of a few hours, we find that scale features such 
as diameters, oblateness and overall shapes are reconstructed with uncertainties 
of a few percent [U [2]. We also quantify the capabilities for reconstructing 
stellar images with non-uniform radiance distributions. We find that hot stars 
(T > 6000° K) containing hot and/or cool locahzed regions (AT - 500° K) as 
small as ~ 0.1 mas can be imaged at short wavelengths (A = 400 nm) [3 . Finally, 
we describe current experimental efforts for measuring intensity correlations in 
a laboratory set-up as well as with a pair of 3 m telescopes separated by 23 m. 
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